Gesture kinetics influence expiratory airflow

Supporting information

Show code setting up the packages

#R-packages
library(papaja)   #for Rmarkdown template for APA(ish) manuscript
library(distill) #for html output
library(ggplot2)    #for plotting
library(gridExtra)  #for plotting in panels
library(knitr)  #for document generation
library(magick)  #for plots
library(cowplot) #for plots
library(plyr) #for revalue
library(tinytex) #for outputting pdfs
library(readr)  #for data reading
library(bookdown) #for crossreferencing
library(rstudioapi) #for setting file paths
library(kableExtra) #for producing nicer tables
library(dplyr) #for data wrangling
library(ggbeeswarm) #for beeswarm plots
library(modelsummary) #for handling model summary tables
library(lme4) #for statistical modeling
library(lmerTest) #for statistical modeling
library(emmeans) #for post-hoc analysis
library(broom.mixed) #for turning fitted models into tidy data frames

#for running python in Rmarkdown
#Sys.setenv("RETICULATE_PYTHON" = "C:/ProgramData/Anaconda3/") # put in your python source
#library(reticulate) # for running python code in Rmarkdown
#if issues with recilate (download the developer version)
 #remotes::install_github("rstudio/reticulate", force = TRUE)

#setwd(normalizePath(".."))
r_refs("references.bib")
#curfol <- normalizePath("..") #set working drive to current

This document contains further information on the supporting methods and results of “Arm movements increase acoustic markers of expiratory flow”. For clarity, it is indicated on the right, what sections in the paper the respective part refers to. Code chunks are hidden by default, but can be made visible by clicking “Show code”.

Supporting methods and results

This is a fully computationally reproducible manuscript written in R Markdown. The full data set is available at the Donders Repository (DSC link).

This refers to the section “Supporting information” of the paper.

Show code loading and preparing the data

#Local data
localfolder <- "C:/Research_Projects/veni_local/exp100/" #THIS NEEDS TO BE CHANGED

#lets set all the data folders
rawd      <- paste0(localfolder, 'Trials/')    #raw trial level data
procd     <- paste0(localfolder, 'Processed/triallevel/') #processed data folder
procdtot  <- paste0(localfolder, 'Processed/complete_datasets/') #processed data folder
daset     <- paste0(localfolder, 'Dataset/')   #dataset folder
meta      <- paste0(localfolder, 'Meta/')   #Meta and trial data

#load in the trialinfo data (containing info like condition order, trial number etc)
trialf    <- list.files(meta, pattern = 'triallist*')
triallist <- data.frame()
for(i in trialf)
{
  tr        <- read.csv(paste0(meta, i), sep = ';')
  tr$ppn    <- parse_number(i)
  triallist <- rbind.data.frame(triallist,tr)
}

#make corrections to trial list due to experimenter errors
  # p10 trial 21 run with a weight
triallist$weight_condition[triallist$ppn==10 & triallist$trial=="21"] <- 'weight' #triallist$trialindex[triallist$ppn==10 & triallist$trial=="21"]
  # p14 double check trial "51", "52", "53" were all weight condition
triallist$weight_condition[triallist$ppn==14 & triallist$trial%in% c("51", "52", "53")] <- 'weight' #triallist$trialindex[triallist$ppn==14 & triallist$trial%in% c("51", "52", "53")]


#load in the metainfo data (containing info like gender, handedness, etc)
mtf     <- list.files(meta, pattern = 'bodymeta*')
mtlist  <- data.frame()
for(i in mtf)
{
  tr      <- read.csv(paste0(meta, i), sep = ';')
  mtlist  <- rbind.data.frame(mtlist,tr)
}
mtlist$ppn <- parse_number(as.character(mtlist$ppn))

Supporting methods

This study has been approved by the Ethics Committee Social Sciences (ECSS) of the Radboud University (reference nr.: 22N.002642).

Show code adding participant and condition information

# add some info about the dataset
participants  <- c(unique(mtlist$ppn))
numpart       <- length(participants)

# number & percentage female & male participants
part_fem <- round(sum(mtlist$sex=='f'), 1)
perc_fem  <- round((sum(mtlist$sex=='f')/numpart)*100, 1)
part_male <- round(sum(mtlist$sex=='m'), 1)
perc_male <- round((sum(mtlist$sex=='m')/numpart)*100, 1)

# body info
m_age           <- round(mean(mtlist$age), 1)
sd_age          <- round(sd(mtlist$age), 1)
m_weight        <- round(mean(mtlist$weight), 1)
sd_weight       <- round(sd(mtlist$weight), 1)
m_height        <- round(mean(mtlist$height), 1)
sd_height       <- round(sd(mtlist$height), 1)
bmi             <- mtlist$weight/((mtlist$height/100)^2)
m_bmi           <- round(mean(bmi), 1)
sd_bmi          <- round(sd(bmi), 1)

# other body measurements
  # Function to calculate the mean of a character vector
  mean_of_char_vector <- function(vec) {
    mean(as.numeric(vec))}

skinfold  <- strsplit(mtlist$upperarmfold, '_')
# Applying the function to each element of the list and storing the results in a new vector
skinfs    <- sapply(skinfold, mean_of_char_vector)
m_skin    <- round(mean(skinfs), 1)
sd_skin   <- round(sd(skinfs), 1)
# we also measure upper/under arm length (but will leave it at this)

# read data
tr_wd <- read.csv(paste0(procdtot, 'fulldata.csv'))
tr_wd <- subset(tr_wd, trialindex > 8 & vocal_condition=='expire') #only real trials and only expiration
numtrials = nrow(tr_wd)
perc_w = round((sum(tr_wd$weight_condition=='weight')/numtrials)*100, 1) #percentage weight
perc_pas = round((sum(tr_wd$movement_condition=='no movement')/numtrials)*100, 1) #percentage passive
perc_int = round((sum(tr_wd$movement_condition=='internal rotation')/numtrials)*100, 1) #percentage internal rotation
perc_ext = round((sum(tr_wd$movement_condition=='external rotation')/numtrials)*100, 1) #percentage external rotation
perc_flex = round((sum(tr_wd$movement_condition=='flexion')/numtrials)*100, 1) #percentage flexion
perc_extens = round((sum(tr_wd$movement_condition=='extension')/numtrials)*100, 1) #percentage extension

Experimental design

This study concerned a two-level wrist-weight manipulation (no weight vs. weight), a two-level within-subject vocalization condition (expire vs. vocalize), and a five-level within-subject movement condition (‘no movement’, ‘extension’, ‘flexion’, ‘external rotation’, ‘internal rotation’). With 4 trial repetitions over the experiment, we yield 80 (2 weights x 2 vocalizations x 5 movements x 4 repetitions) trials per participant. Trials were blocked by weight condition and vocalization condition (so that weights and task did not switch from trial to trial). Within blocks all movement conditions were randomized.

Participants

For the current pre-registered confirmatory experiment, as planned and supported by a power analysis (see preregistration), we collected N = (17) participants: 7 female, 10 male, M (SD) age = 28.50 (6.50), M (SD) body weight = 72.10 kg (10.20), M (SD) body height = 175.10 cm (8.50), M (SD) BMI = 23.40 (2.20), M (SD) triceps skinfold = 19.10 mm (4.30).

Exclusions and deviations from pre-registration

We also performed the experiment with one other participant, but due to an issue with LSL streaming, this dataset could not be synchronized and was lost. Furthermore, due to running over time one participant had to terminate the study earlier about half way through. Note further, that in our pre-registration we wanted to admit participants with a BMI lower than 25, but since participants were difficult to recruit we accepted three participants with slightly higher BMIs too (max BMI of the current dataset: 27.10). Participants were all able-bodied and did not have any constraints in performing the task. Finally, as stated in the pre-registration we will only report results on the expiration trials (and not the vocalization-only trials). This means that for this report in total we have N = 17 participants, with 636 analyzable trials, with balanced conditions: weight = 50.50%, no movement = 20.00%, internal rotation = 20.00 %, external rotation = 20.00%, flexion = 20.00 %, extension = 20.00 %.

Measurements and equipment

Body measurements

To enable future analyses of possible factors modulating individual-specific body properties, we collect some basic information about body properties. Namely, weight, under arm length, upper arm length, triceps skinfold, and upper arm circumference.

Experiment protocol

The experiment was coded in Python using functions from PsychoPy. The experiment was controlled via a Brainvision Button Box (Brain Products GmbH, Munich, Germany), which was also streaming its output to the data collection PC unit.

Wrist weights

To manipulate the mass set in motion, we apply a wrist weight. We use a TurnTuri sports wrist weight of 1 kg.

Video and kinematics

The participants are recorded via a videocamera (Logitech StreamCam), sampling at 60 frames per second. We used Mediapipe (Lugaresi et al. 2019) to track the skeleton and facial movements, which is implemented in Masked-piper which we also use for masking the videos (Owoyele et al. 2022). The motion-tracked skeleton, specifically the wrist of the dominant hand, is used to estimate movement initiation, peak speed, and the end of the movement. The motion tracking is, however, only used for determining movement windows, and is not of central concern.

Muscle activity (Surface ElectroMyography: sEMG)

We measured sEMG using a wired BrainAmp ExG system (Brain Products GmbH, Munich, Germany). Disposable surface electrodes (Kendall 24mm Arbo H124SG) were used, and for each of the four muscle targets we had 3 (active, reference, ground) electrodes (12 electrodes total). The sEMG system sampled at 2,500 Hz (for post-processing filters see below).
For an overview of the electrode attachments, see Fig. S1. We prepare the skin surface for EMG application with a scrub gel (NuPrep) followed by cotton ball swipe with alcohol (Podior 70 %). Active and reference electrodes were attached with a 15mm distance center to center.

Figure 1: Overview sEMG target muscles. Active (a) and reference (r), and ground (g) sEMG electrodes, for each muscle target.

We attached electrodes for focal muscles which directly participate in the internal (pectoralis major) and external rotation (infraspinatus) of the humerus. Electrodes were applied for focal muscles ipsilaterally (relative to the dominant hand). We attached electrodes to the muscle belly of the clavicular head of the pectoralis major, with a ground electrode on the clavicle on the opposite side.
We also attached electrodes for postural muscles which will likely anticipate and react to postural perturbations due to upper limb movements. Since these muscles should act in the opposite direction of the postural perturbation of the dominant hand, we applied electrodes contralaterally to the dominant hand. We attach electrodes to the rectus abdominis, with a ground electrode on the iliac crest on the opposite side. We also attached electrodes to the erector spinae muscle group (specifically, the iliocostalis lumborum).

Ground reaction measurements

We used an inhouse-built 1m² balance board with vertical pressure sensors. The sensors were derived and remodified from the Wii-Balance board sensors. The sampling rate was 400 Hz. The system was time-locked within millisecond accuracy, and has a spatial accuracy of several sub-millimeters. A national instruments card, USB-62221 performed the A/D conversion and was connected via USB to the PC.

Acoustics

To ensure proper acoustic intensity measurements we used a headset microphone; MicroMic C520 (AKG, Inc.) headset condenser cardioid microphone sampling at 16 kHz. The gain levels of the condenser power source were set by the hardware (and could not be changed).

Here are five audio examples from a male participant producing the exhalation in different conditions without added weight. You might have to set your volume a bit higher than usually for them.

Show code setting up the audio examples

library(htmltools)

audioexamples <- paste0(basefolder, "/AudioExamples/")

audioexample_nomovement       <- paste0(audioexamples, "expire_sample1_nomovement.wav")
audioexample_internalrotation <- paste0(audioexamples, "expire_sample2_internalrotation.wav")
audioexample_flexion          <- paste0(audioexamples, "expire_sample3_flexion.wav")
audioexample_externalrotation <- paste0(audioexamples, "expire_sample4_externalrotation.wav")
audioexample_extension        <- paste0(audioexamples, "expire_sample5_extension.wav")

html_tag_audio <- function(file, type = c("wav")) {
  type <- match.arg(type)
  htmltools::tags$audio(
    controls = NA,
    htmltools::tags$source(
      src = file,
      type = glue::glue("audio/{type}", type = type)
    )
  )
}

The first example is from the no movement condition:

The second example is from the internal rotation condition:

The third example is from the external rotation condition:

The fourth example is from the extension condition:

The fifth example is from the flexion condition:

Recording setup and synchronization

We use LabStreamLayer that provides a uniform interface for streaming different signals along a network, where a common time stamp for each signal ensures sub-millisecond synchronization. We used a Linux system to record and stream the microphone recordings. Additionally a second PC collected video, and streamed ground reaction forces, and EMG. A data collection PC collected the audio, ground reaction force, and EMG streams and stored the output in XDF format for efficient storing of multiple time-varying signals. The video data was synchronized by streaming the frame number to the LSL recorder, allowing us to match up individual frames with the other signals (even when a frame is dropped).

Procedure

Participants are admitted to the study based on exclusion criteria and sign an informed consent. We ask participants to take off their shoes and we proceed with the body measurements, while instructing the participant about the nature of the study. After body measurements, we apply the surface EMG. We prepared the muscle site with rubbing gel and alcohol, and active/reference electrodes were placed with a distance of 15 mm from each other’s center. See Fig. S1 for the sEMG electrode locations. The procedures up to the start of the experiment take about 20 minutes or less in total.

Upon the start of the experiment participants take a standing position on the force platform. The experiment commences with calibration and practice trials. First, 10 seconds of silent breathing with body movements are recorded. Then participants are asked to take a maximum inspiration followed by a maximum expiration to measure signal conditions under respiratory boundary conditions. Then, for the practice trials, each movement was practiced with expiring and vocalization while performing the movement conditions, and the participant is introduced to wearing the wrist weight of 1 kg. After practice trials, participants performed 80 blocked trials.

For each (practice) trial participants are closely guided by the information on the monitor. Firstly, participants are shown the movement to be performed for the trial, and have to prompt the experimenter that they are ready to continue. Then participants are instructed to adopt the start position of the movement, which is a 90 degree elbow flexion, with either an externally rotated humerus (start position for internal rotation), or a non-rotated humerus with the wrist in front of the body (rest position for the other movement conditions). For the no movement condition participants are asked to rest their arms alongside their body. Upon trial start, participants are asked to inhale deeply with a timer counting down from 4 seconds. Then, participants are asked to ‘vocalize’ or ‘expire’, with a screen appearing after 3 seconds to perform the movement with visual guidance to where the movement end position is so that participants are reminded of the movement. After an additional 4 seconds the trial ends, which allows more than enough time to perform the movement and stabilize vocalization/expiration after the perturbation. In the no movement condition, a prompt is given to maintain one’s posture.

Preprocessing of the data streams

EMG

To reduce heart rate artifacts we apply a common method (Drake and Callaghan 2006) of high-pass filtering the signal at 30 Hz using a zero-phase 4th order Butterworth filter. We then full-wave rectified the EMG signal and applied a zero-phase low-pass 4th order Butterworth filter at 20 Hz. When filtering any signal we pad the signals to avoid edge effects. We normalized the EMG signals within participants before submitting to analyses.

Show code preprocessing EMG data

# load in the data and save it
library(signal)
library(ggrepel)

#####################preprocessing functions
butter.it <- function(x, samplingrate, order, cutoff, type = "low")
{
  nyquist <- samplingrate/2
  xpad <- c(rep(0, 1000), x, rep(0,1000)) #add some padding
  bf <- butter(order,cutoff/nyquist, type=type) 
  xpad <- as.numeric(signal::filtfilt(bf, xpad))
  x <- xpad[1000:(1000+length(x)-1)] #remove the padding
}

#rectify and high and low pass EMG signals
reclow.it <- function(emgsignal, samplingrate, cutoffhigh, cutofflow)
  {
    #high pass filter
    output <- butter.it(emgsignal, samplingrate=samplingrate, order=4, cutoff = cutoffhigh,
                        type = "high")
    #rectify and low pass
    output <<- butter.it(abs(output), samplingrate=samplingrate,  order=4, cutoff= cutofflow,
                        type = "low")
}
###########################################

####################################LOAD IN EMG DATA
EMGdata <- list.files(rawd, pattern = '*Dev_1.csv')
####################################SET common colors for muscles
col_pectoralis    <- '#e7298a'
col_infraspinatus <- '#7570b3'
col_rectus        <- '#d95f02'
col_erector       <- '#1b9e77'
colors_mus <- c("pectoralis major" = col_pectoralis, "infraspinatus" = col_infraspinatus, "rectus abdominis" = col_rectus, "erector spinae" = col_erector)

#show an example
emgd <- read.csv(paste0(rawd, EMGdata[822]))
colnames(emgd) <- c('time_s', 'infraspinatus', 'rectus_abdominis', 'pectoralis_major', 'erector_spinae', 'empty1', 'empty2', 'empty3', 'empty4', 'respiration_belt', 'button_box')
samplingrate <- 1/mean(diff(emgd$time_s))
nyquistemg <- samplingrate/2

#example without strong high pass filter
emgd$time_s <- emgd$time_s-min(emgd$time_s)
emgd$pectoralis_major_sm <- reclow.it(emgd$pectoralis_major, samplingrate= samplingrate, cutoffhigh = 2, cutofflow = 20)
emgd$infraspinatus_sm <- reclow.it(emgd$infraspinatus, samplingrate= samplingrate, cutoffhigh = 2, cutofflow = 20)
emgd$erector_spinae_sm <- reclow.it(emgd$erector_spinae, samplingrate= samplingrate, cutoffhigh = 2, cutofflow = 20)
emgd$rectus_abdominis_sm <- reclow.it(emgd$rectus_abdominis, samplingrate= samplingrate, cutoffhigh = 2, cutofflow = 20)

#example trial
unfilteredEMG <- ggplot(emgd[seq(1, nrow(emgd),by=5),], aes(x=time_s))+geom_path(aes(y=infraspinatus_sm, color = 'infraspinatus'))+
  geom_path(aes(y=pectoralis_major_sm, color = 'pectoralis major'))+
  geom_path(aes(y=rectus_abdominis_sm, color = 'rectus abdominis'))+
  geom_path(aes(y=erector_spinae_sm, color = 'erector spinae'))+ 
  geom_text(aes(x = 0.5, y = 600), color = "black", label = "Heart rate peak 1", angle=45, size = 2) +
  geom_text(aes(x = 1.4, y = 600), color = "black", label = "Heart rate peak 2", angle=45, size = 2) +
  geom_text(aes(x = 4.4, y = 600), color = "black", label = "Heart rate peak 5", angle=45, size = 2) +
  scale_color_manual(values = colors_mus)
unfilteredEMG <- unfilteredEMG +
  labs(x = "time in seconds", y = "EMG rectified") +
  xlab('time in seconds') +
  ylab('EMG rectified') +
  ylim(0, 800) +
  theme_cowplot(12) +
  theme(legend.position="none")


#example without strong high pass filter
emgd$pectoralis_major_sm_h <- reclow.it(emgd$pectoralis_major, samplingrate= samplingrate, cutoffhigh = 30, cutofflow = 20)  
emgd$infraspinatus_sm_h  <- reclow.it(emgd$infraspinatus, samplingrate= samplingrate, cutoffhigh = 30, cutofflow = 20) 
emgd$erector_spinae_sm_h  <- reclow.it(emgd$erector_spinae, samplingrate= samplingrate, cutoffhigh = 30, cutofflow = 20) 
emgd$rectus_abdominis_sm_h  <- reclow.it(emgd$rectus_abdominis, samplingrate= samplingrate, cutoffhigh = 30, cutofflow = 20) 

#example trial
filteredEMG <- ggplot(emgd[seq(1, nrow(emgd),by=5),], aes(x=time_s)) +
  geom_path(aes(y=infraspinatus_sm_h, color = 'infraspinatus')) +
  geom_path(aes(y=pectoralis_major_sm_h, color = 'pectoralis major')) +
  geom_path(aes(y=rectus_abdominis_sm_h, color = 'rectus abdominis')) +
  geom_path(aes(y=erector_spinae_sm_h, color = 'erector spinae')) +
  scale_color_manual(values = colors_mus)
filteredEMG <- filteredEMG +
  labs(x = "time in seconds", y = "EMG rectified", color = "Legend") +
  xlab('time in seconds') +
  ylab('EMG rectified') +
  theme_cowplot(12)+
  theme(
    legend.position = c(1, 1),
    legend.justification = c("right", "top"),
    legend.box.just = "right",
    legend.margin = margin(4, 4, 4, 4)
    )

Figure 2: Example of smoothing settings for EMG signals. The upper panel shows the raw rectified high-pass filtered EMG, and the lower panel shows the low-pass filtered data to reduce artifacts of heart rate. This example shows an internal rotation trial, where we successfully retrieve the peak in pectoralis major that internally rotates the arm.

Ground reaction forces

We upsampled the balanceboard from 400 Hz to 2,500 Hz. We then applied a zero-phase low-pass 20 Hz 2nd order Butterworth filter to the padded signals. As a key measure for postural perturbation we computed the change in 2D magnitude (L2 norm of the center of pressure x and y) in center of pressure (hereafter COPc).

Acoustics

For acoustics we extract the smoothed amplitude envelope (hereafter envelope). For the envelope we apply a Hilbert transform to the waveform signal, then take the complex modulus to create a 1D time series, which is then resampled at 2,500 Hz, and smoothed with a Hann filter based on a Hanning Window of 12 Hz. We normalize the amplitude envelope signals within participants before submitting to analyses.

Show code extracting acoustics

library(rPraat) #for reading in sounds
library(dplR)   #for applying Hanning
library(seewave) #for signal processing general
library(wrassp) #for acoustic processing

##################### MAIN FUNCTION TO EXTRACT SMOOTHED ENVELOPE ###############################
amplitude_envelope.extract <- function(locationsound, smoothingHz, resampledHz)
{
  #read the sound file into R
  snd <- rPraat::snd.read(locationsound)
  #apply the hilbert on the signal
  hilb <- seewave::hilbert(snd$sig, f = snd$fs, fftw =FALSE)
  #apply complex modulus
  env <- as.vector(abs(hilb))
  #smooth with a hanning window
  env_smoothed <- dplR::hanning(x= env, n = snd$fs/smoothingHz)
  #set undeterminable at beginning and end NA's to 0
  env_smoothed[is.na(env_smoothed)] <- 0
  #resample settings at desired sampling rate
  f <- approxfun(1:(snd$duration*snd$fs),env_smoothed)
  #resample apply
  downsampled <- f(seq(from=0,to=snd$duration*snd$fs,by=snd$fs/resampledHz))
  #let function return the downsampled smoothed amplitude envelope
  return(downsampled[!is.na(downsampled)])
}

Data aggregation

All signals were sampled at, or upsampled to, 2,500 Hz. Then we aggregated the data by aligning time series in a combined by-trial format to increase ease of combined analyses. We linearly interpolated signals when sample times did not align perfectly. For the amplitude envelope, we wanted to transform it to dB values. As a reference value, we wanted to use the lab under quiet conditions. So we selected participant 1’s practice trial 1 while not producing any noise for that (from second 2 to second 4) and extracted the mean power from that.

Show code aggregating the data

library(pracma)
library(zoo)

overwrite = T

#lets loop through the triallist and extract the files
triallist$uniquetrial <- paste0(triallist$trialindex, '_', triallist$ppn) 
  #remove duplicates so that we ignore repeated trials
howmany_repeated <- sum(duplicated(triallist$uniquetrial))
triallist <- triallist[!duplicated(triallist$uniquetrial),]

#get reference value first, so we can turn env into envdb
# referenceenv <- amplitude_envelope.extract(locationsound = paste0(localfolder, "Trials/p1_trial_1_Micdenoised.wav"), smoothingHz =  12, resampledHz = 2500)
# referenceenvcut <- referenceenv[5000:10000]
# referenceenvmean <- mean(referenceenvcut)

for(ID in triallist$uniquetrial)
{
  #print(paste0('working on ', ID))
  #debugging
  #ID <- triallist$uniquetrial[27] #I commented that out bc it was so many lines
  
  triallist_s <- triallist[triallist$uniquetrial==ID,]
  ##load in row triallist
  pp <- triallist_s$ppn
  pps <- triallist_s$ppn
  if((pp == '41') | (pp == '42')){pps <- '4'}
  gender <- mtlist$sex[mtlist$ppn==pps]
  hand   <- mtlist$handedness[mtlist$ppn==pps]
  triali <- triallist_s$trialindex
  
  if(!file.exists(paste0(procd, 'pp_', pp, '_trial_', triali, '_aligned.csv')) | overwrite == TRUE)
  {
  ######load in EMG and resp belt
    EMGr <- read.csv(paste0(rawd, 'p', pp, '_trial_',triali, '_BrainAmpSeries-Dev_1.csv'))
    channels <- c('time_s', 'infraspinatus', 'rectus_abdominis', 'pectoralis_major', 'erector_spinae', 'empty1', 'empty2', 'empty3', 'empty4', 'respiration_belt', 'button_box')
    colnames(EMGr) <- channels
    #only keep the relevant vectors
    EMGr <- cbind.data.frame(EMGr$time_s,EMGr$respiration_belt, EMGr$infraspinatus,
                             EMGr$rectus_abdominis,EMGr$pectoralis_major, EMGr$erector_spinae)
    colnames(EMGr) <- channels[c(1, 10, 2, 3, 4, 5)] #rename
    EMGsamp <- 1/mean(diff(EMGr$time_s-min(EMGr$time_s)))
    #smooth
      #note the respiration belt is not installed now, but if we resolve technical issues we might add it
    EMGr$respiration_belt <- butter.it(EMGr$respiration_belt, samplingrate = EMGsamp, order = 2, cutoff = 30)
      #smooth EMG with described settings
    EMGr[,3:6] <- apply(EMGr[,3:6], 2, function(x) reclow.it(scale(x, scale = FALSE, center =TRUE), samplingrate= EMGsamp, 
                                                             cutoffhigh = 30, cutofflow = 20)) #CHANGECONF we center before we smooth now
  #########load in acoustics
  locsound <- paste0(rawd, 'p', pp, '_trial_',triali, '_Micdenoised.wav')
  if(!file.exists(locsound)){print(paste0('this mic file does not exist: ', locsound))}
  if(file.exists(locsound)){
    
  
  env <- amplitude_envelope.extract(locationsound=locsound, smoothingHz =  12, resampledHz = 2500)
  #envdb <- 10 * log10(env/referenceenvmean)
  acoustics <- cbind.data.frame(env) #add envdb here if you want it
  acoustics$f0 <- NA
  acoustics$time_s <- seq(0, (nrow(acoustics)-1)*1/2500, 1/2500)
  duration_time <- (max(acoustics$time_s)-min(acoustics$time_s)) #note down duration time of this trial
  #########load in motion
  mt <- read.csv(paste0(rawd, 'p', pp, '_trial_',triali, '_video_raw_body.csv'))
  if(hand == 'r'){mov <- cbind.data.frame(mt$X_RIGHT_WRIST, mt$Y_RIGHT_WRIST, mt$Z_RIGHT_WRIST)}  #CHANGECONF change we also take depth
  if(hand == 'l'){mov <- cbind.data.frame(mt$X_LEFT_WRIST, mt$Y_LEFT_WRIST, mt$Z_RIGHT_WRIST)}    #CHANGECONF we also take depth
  colnames(mov) <- c('x_wrist', 'y_wrist', 'z_wrist')
  mov$y_wrist <- mov$y_wrist*-1 #flip axes so that higher up means higher values (check)
  #filter movements
  mov[,1:3] <- apply(mov[,1:3], 2, FUN=function(x) butter.it(x, order=4, samplingrate = 60, cutoff = 10, type='low'))
  
  #########load in balanceboard
  bb <- read.csv(paste0(rawd, 'p', pp, '_trial_',triali, '_BalanceBoard_stream.csv'))
  colnames(bb) <- c('time_s', 'left_back', 'right_forward', 'right_back', 'left_forward')
  bbsamp <- 1/mean(diff(bb$time_s - min(bb$time_s)))
  bb[,2:5] <- apply(bb[,2:5], 2, function(x) sgolayfilt(x, 5, 51))
  COPX  <- (bb$right_forward+bb$right_back)-(bb$left_forward+bb$left_back)
  COPY  <- (bb$right_forward+bb$left_forward)-(bb$left_back+bb$left_back)
  bb$COPXc <- sgolayfilt(c(0, diff(COPX)), 5, 51)
  bb$COPYc <-  sgolayfilt(c(0, diff(COPY)), 5, 51)
  bb$COPc <- sqrt(bb$COPXc^2 + bb$COPYc^2)
  #########combine everything
    #combined acoustics and EMG(+resp)
  EMGr$time_ss <- EMGr$time_s-min(EMGr$time_s)
  ac_emg <- merge(acoustics, EMGr, by.x = 'time_s', by.y = 'time_ss', all=TRUE)
  ac_emg[,4:9] <- apply(ac_emg[,4:9], 2, function(y) na.approx(y, x=ac_emg$time_s, na.rm = FALSE))
  ac_emg <- ac_emg[!is.na(ac_emg$env),]
  ac_emg$time_s <- ac_emg$time_s + min(ac_emg[,4],na.rm=TRUE) #get the original time
  ac_emg <- ac_emg[!is.na(ac_emg$time_s),-4] #remove the centered time variable, and remove time NA at the end
    #now add balance board
  ac_emg_bb <- merge(ac_emg, bb[,c(1, 6:8)], by.x='time_s', by.y='time_s', all= TRUE)
  ac_emg_bb[,9:11] <- apply(ac_emg_bb[,9:11], 2, function(y) na.approx(y, x=ac_emg_bb$time_s, na.rm = FALSE))
  ac_emg_bb <- ac_emg_bb[!is.na(ac_emg_bb$env), ]
  ac_emg_bb <- ac_emg_bb[!is.na(ac_emg_bb$time_s),-4]
    #now add the final motion (which is thereby upsampled to 2500 hertz)
  start_time <- min(ac_emg_bb$time_s,na.rm=TRUE)
  mov$time_s <-   start_time+seq(0, duration_time-(duration_time /nrow(mov)), by= duration_time /nrow(mov))
  ac_emg_bb_m <- merge(ac_emg_bb, mov, by.x='time_s', by.y='time_s', all = TRUE)
  ac_emg_bb_m[,11:13] <- apply(ac_emg_bb_m[,11:13], 2, function(y) na.approx(y, x=ac_emg_bb_m$time_s, na.rm = FALSE))
  ac_emg_bb_m <- ac_emg_bb_m[!is.na(ac_emg_bb_m$env), ]
  #now write the files away to a processed folder
  write.csv(ac_emg_bb_m, paste0(procd, 'pp_', pp, '_trial_', triali, '_aligned.csv'))
    }
  }
}

Data sharing & Privacy

Video data is deidentified using the masked-piper tool to mask faces and body while maintaining kinematic information (Owoyele et al. 2022).

Overview data analyses

Note that there may be a general decline in the amplitude of the exhalation during a trial (see Fig. S3, top panel). Unlike in the vocalization condition, where a decline is to be expected, as the subglottal pressure falls when the lungs deflate, this is not always the case here. Strategies for producing a sustained exhalation differ and as a consequence, the amplitude may sometimes decrease, increase, or even remain stable throughout. To cater for that and to quantify deviations from stable exhalations, we therefore detrend the amplitude envelope time series, so as to assess positive or negative peaks relative to this trend line. For the envelope, muscle activity, and the change in center of pressure we will measure the global maxima happening within the analyses window (i.e., within a trial we take a local maximum occurring between movement onset and offset). We will analyze positive and negative peaks within the movement window separately.

Show code scaling data within participant

#We have to do some within participant scaling, such that
  #all EMG's are rescaled within participants
fs <-  list.files(procd)
fs <- fs[!fs%in%list.files(procd, pattern = 'zscal_*')]

#set to true if you want to overwrite files (to check [modified] code)
overwrite = T


if((length(list.files(procd, pattern = 'zscal_*')) == 0) | (overwrite==TRUE))
{
  mr <- data.frame()
  for(f in fs)
  {
    temp <- read.csv(paste0(procd, f))
    temp$pp <- parse_number(f)
    temp$tr <- f
    mr <- rbind.data.frame(mr, temp)
  }
  
  #set ppn to 4
  temp$pp[temp$pp==41 | temp$pp==42] <- 4
  
  #center the muscle measurements for each trial (due to drift)
  mr$infraspinatus <- ave(mr$infraspinatus, mr$tr, FUN = function(x){x<-x-mean(x, na.rm=TRUE)})       
  mr$rectus_abdominis <- ave(mr$rectus_abdominis, mr$tr, FUN = function(x){x<-x-mean(x, na.rm=TRUE)}) 
  mr$pectoralis_major <- ave(mr$pectoralis_major, mr$tr, FUN = function(x){x<-x-mean(x, na.rm=TRUE)}) 
  mr$erector_spinae <- ave(mr$erector_spinae, mr$tr, FUN = function(x){x<-x-mean(x, na.rm=TRUE)})     

  
  #normalize signals by participant
  mr$infraspinatus <- ave(mr$infraspinatus, mr$pp, FUN = function(x)scale(x, center = FALSE))       
  mr$rectus_abdominis <- ave(mr$rectus_abdominis, mr$pp, FUN = function(x)scale(x, center = FALSE)) 
  mr$pectoralis_major <- ave(mr$pectoralis_major, mr$pp, FUN = function(x)scale(x, center = FALSE)) 
  mr$erector_spinae <- ave(mr$erector_spinae, mr$pp, FUN = function(x)scale(x, center = FALSE))     
  mr$env_z <- ave(mr$env, mr$pp, FUN = function(x)scale(x, center = FALSE))
# we cannot z-scale envdb here, because there are -Inf (resulting from 0 in env), as a consequence, we will do this in the next step
  
  #save the scaled variables in new objects
  for(f in unique(mr$tr))
  {
    sb <- mr[mr$tr==f, ] #subset
    #center motion per trial
    sb$x_wrist <- sb$x_wrist-mean(sb$x_wrist,na.rm=TRUE)
    sb$y_wrist <- sb$y_wrist-mean(sb$y_wrist,na.rm=TRUE)
    sp <- c(0, sqrt(diff(sb$x_wrist)^2+diff(sb$y_wrist)^2))
    sp[is.na(sp)] <- 0 #there are some trailing NA's which we fill with 0 
    sb$speed <- butter.it(sp, samplingrate = 2500, order = 4, cutoff = 5)
    sb$time_s   <- sb$time_s #this is the time
    sb$time_s_c <- sb$time_s-min(sb$time_s)# accumulating at trial start (resyncing) resynced with the psychopy software
    write.csv(sb, paste0(procd, 'zscal_', f))
  }
}

Show code producing the plot

ptest=8
ex2 <- triallist[triallist$vocal_condition=='expire' & triallist$movement_condition=='extension_stop'& triallist$ppn ==ptest, ]
ex2 <- ex2[1, ]
#ex2$trialindex
exts2 <- read.csv(paste0(procd, 'zscal_pp_', ex2$ppn, '_trial_', ex2$trialindex, '_aligned.csv'))
#exts2 <- dplyr::rename(exts2, rectus_abdominis = rectus_abdominus)
  exts2 <- subset(exts2, time_s_c > 1 & time_s_c <5)
  exts2 <- exts2[seq(1, nrow(exts2), by = 5), ]    #downsample by a 5th
  exts2$time_ms <- round(exts2$time_s_c*1000)
  fp <- pracma::findpeaks(exts2$speed)
  fp <- fp[which.max(fp[,1]),]
  startmov <- exts2$time_ms[fp[3]]
  endmov   <- exts2$time_ms[fp[4]]
  ps       <- exts2$time_ms[fp[2]]
#apply movement rule to centralize trials
#here all processing steps are performed, and data aligned
exampleenvelope <- ggplot(exts2, aes(x=time_ms)) +
  geom_path(aes(y=env), size =1) +
  theme_cowplot(12) + 
  ggtitle('Example trial extension') +
  ylab('norm. env') +
  geom_vline(xintercept=startmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=endmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=startmov-500, color = 'black', size= 2, alpha=.3) +
  geom_smooth(aes(y=env), method='lm', linetype = 'dashed',color = 'black') +
  geom_vline(xintercept=ps, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept= 3000) +
  xlim(1000, 5000) +
  #ylim(27.5, 35) + #only needed for envdb
  xlab('time (ms)')

examplemovement <- ggplot(exts2, aes(x=time_ms)) +
  geom_path(aes(y=round(speed*2500*100)), color = 'black') +
  theme_cowplot(12) +
  geom_vline(xintercept=startmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=endmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=startmov-500, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept=ps, color = 'black', size= 2, alpha=.3) +
  ylab("cm per second") +
  geom_vline(xintercept=3000, color = 'black', size= 1) +
  geom_vline(xintercept= 3000) +
  xlim(1000, 5000)+xlab('time (ms)')

exampleCOP <- ggplot(exts2, aes(x=time_ms)) +
  geom_path(aes(y=COPc)) +
  theme_cowplot(12) +
  geom_vline(xintercept=startmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=endmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=startmov-500, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept=ps, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept= 3000) +xlim(1000, 5000) +
  xlab('time (ms)')

exampleEMG <- ggplot(exts2, aes(x=time_ms)) +
  geom_path(aes(y=infraspinatus, color = 'infraspinatus')) +
  geom_path(aes(y=pectoralis_major, color = 'pectoralis major')) +
  geom_path(aes(y=rectus_abdominis, color = 'rectus abdominis')) +
  geom_path(aes(y=erector_spinae, color = 'erector spinae')) +
  geom_vline(xintercept=startmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=endmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=startmov-500, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept=ps, color = 'black', size= 2, alpha=.3) +
  scale_color_manual(values = colors_mus) +
  labs(y="sEMG") +
  theme_cowplot(12) +
  theme(legend.position = c(0.8, 0.6)) +
  geom_vline(xintercept= 3000) +
  xlim(1000, 5000) +
  xlab('time (ms)')

Figure 3: Combined time series. One example trial and the associated signals are shown, for the internal rotation movement condition. At time = 0 s, the prompt is given to the participant to expire. We determine a detrending line using linear regression for the 1 to 5 s after the expiration prompt. Note, at 3 s (3000 ms), there is a movement prompt. However, we determine our window where we assess peaks in signals at 500 ms before and after the movement onset/offset (using peakfinding function on the 2D speed time series of the wrist). In these trials, the analysis window is given in gray dashed bars, which is 500 ms after and before movement onset/offset.

Show code setting up the data frame

overwrite = T #change this back to FALSE

if(overwrite)
{
tr_wd <- triallist
tr_wd$min_amp_c_around_move <- NA
tr_wd$min_amp_time_around_move <- NA
tr_wd$max_amp_c_around_move <- NA
tr_wd$max_amp_time_around_move  <- NA
tr_wd$max_pectoral <- NA
tr_wd$max_infra<- NA
tr_wd$max_rectus<- NA
tr_wd$max_erector<- NA
tr_wd$max_COPxc <- NA
tr_wd$max_COPyc <- NA
tr_wd$max_COPc  <- NA

#lets loop through the triallist and extract the files
tr_wd$uniquetrial <- paste0(tr_wd$trialindex, '_', triallist$ppn) 
  #remove duplicates so that we ignore repeated trials #NOTE in total we registered 25 repeats
howmany_repeated <- sum(duplicated(tr_wd$uniquetrial))
tr_wd <- tr_wd[!duplicated(tr_wd$uniquetrial),]

#make one large dataset
tsl <- data.frame()

for(ID in tr_wd$uniquetrial)
{
   #debugging
  #ID <- tr_wd$uniquetrial[27]
  tr_wd_s <- tr_wd[tr_wd$uniquetrial==ID,]
  
  ##load in row tr_wd
  pp <- tr_wd_s$ppn
  gender <- mtlist$sex[mtlist$ppn==pp]
  hand   <- mtlist$handedness[mtlist$ppn==pp]
  triali <- tr_wd_s$trialindex
  if(file.exists(paste0(procd, 'zscal_pp_', pp, '_trial_', triali, '_aligned.csv')))
  {
  ts <- read.csv(paste0(procd, 'zscal_pp_', pp, '_trial_', triali, '_aligned.csv'))
  ts$time_ms <- round(ts$time_s*1000)
  ts$time_ms <- ts$time_ms-min(ts$time_ms)
                                              #take the first 1 to 5 seconds
  ts <- ts[(ts$time_ms > 1000) & (ts$time_ms < 5000), ]
  ts$env_z <- pracma::detrend(ts$env_z, tt = 'linear') # detrend the envelope
  #z-scaling envdb now happens here
  # ts <- ts %>%
  #   #dplyr::filter(env != 0) %>%
  #   group_by(pp) %>%
  #   mutate(envdb_z = scale(envdb, center = FALSE)) %>%
  #   ungroup
  # 
  # ts$envdb_z <- pracma::detrend(ts$envdb_z, tt = 'linear') # detrend the envelope
  #movement onset + - 500 milliseconds
  fp <- pracma::findpeaks(ts$speed)
  fp <- fp[which.max(fp[,1]),]
  startmov <- ts$time_ms[fp[3]]-500
  endmov   <- ts$time_ms[fp[4]]+500
  ts$movstart <- ts$time_ms-ts$time_ms[fp[3]] #get the time when the movement starts (which is the beginning of the rise to peak speed)
  
  subts <- ts[(ts$time_ms>startmov) & (ts$time_ms<endmov),] #subset
  subts$peaks <- NA
  indexspeed <- which(subts$time_ms==ts$time_ms[fp[2]])[1] #there are might be multiple adjacent values that have the peak speed
  subts$peaks[indexspeed] <- 'peak speed'
  indexend <- which(subts$time_ms==ts$time_ms[fp[4]])[1] #there are might be multiple adjacent values that have the peak speed
  subts$peaks[indexend] <- 'movement end'
  subts$movement_condition <- tr_wd_s$movement_condition
  subts$vocal_condition <- tr_wd_s$vocal_condition
  subts$weight_condition <- tr_wd_s$weight_condition
  subts$trial_number <-  tr_wd_s$trialindex
  
  
  #bind into big dataset
  tsl <- rbind(tsl, subts)

  #peakspeed
  tp <- subts$time_ms[which(subts$peaks=='peak speed')] #take everything before the peak speed so that we get movement onset

  # #get the subts
  tr_wd$max_amp_c_around_move[tr_wd$uniquetrial==ID] <- max(subts$env_z[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  tr_wd$min_amp_c_around_move[tr_wd$uniquetrial==ID] <- min(subts$env_z[subts$time_ms<tp], na.rm=TRUE) #before the peak speed

  # #get the subts
  # tr_wd$max_amp_c_around_move[tr_wd$uniquetrial==ID] <- max(subts$envdb_z[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  # tr_wd$min_amp_c_around_move[tr_wd$uniquetrial==ID] <- min(subts$envdb_z[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  # 
  # # #also do it for un-zscaled envelope
  # tr_wd$max_amp_c_around_move_unscaled[tr_wd$uniquetrial==ID] <- max(subts$envdb[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  # tr_wd$min_amp_c_around_move_unscaled[tr_wd$uniquetrial==ID] <- min(subts$envdb[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  
  # also collect the times of the peaks
  envcc <- subts$env_z[subts$time_ms<tp]
  timecc <- subts$time_ms[subts$time_ms<tp]
  tr_wd$max_amp_time_around_move[tr_wd$uniquetrial==ID] <- timecc[which.max(envcc)]-tp # time pos peak
  tr_wd$min_amp_time_around_move[tr_wd$uniquetrial==ID] <- timecc[which.min(envcc)]-tp #time negative peak relative to peak speed
  # collect the peaks
  tr_wd$max_pectoral[tr_wd$uniquetrial==ID] <- max(subts$pectoralis_major[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  tr_wd$max_infra[tr_wd$uniquetrial==ID] <- max(subts$infraspinatus[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  tr_wd$max_rectus[tr_wd$uniquetrial==ID] <- max(subts$rectus_abdominis[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  tr_wd$max_erector[tr_wd$uniquetrial==ID] <- max(subts$erector_spinae[subts$time_ms<tp], na.rm=TRUE)#before the peak speed
  tr_wd$max_COPxc[tr_wd$uniquetrial==ID] <- max(subts$COPXc[subts$time_ms<tp], na.rm=TRUE)#before the peak speed
  tr_wd$max_COPyc[tr_wd$uniquetrial==ID] <- max(subts$COPYc[subts$time_ms<tp], na.rm=TRUE)#before the peak speed
  tr_wd$max_COPc[tr_wd$uniquetrial==ID] <- max(subts$COPc[subts$time_ms<tp], na.rm=TRUE)#before the peak speed
  }
}

#write a large dateset to 
#write.csv(tsl, )

#order levels
tr_wd$movement_condition <- recode(tr_wd$movement_condition, no_movement  = "no movement",
                                           extension_stop = "extension",flexion_stop  = "flexion",
                                           external_rotation_stop = "external rotation",
                                           internal_rotation_stop = "internal rotation")
tr_wd$movement_condition <- factor(tr_wd$movement_condition, 
                                   levels = c("no movement", "extension", "flexion", 
                                              "external rotation", "internal rotation"))

#maybe write to file?
write.csv(tr_wd, paste0(procdtot, 'fulldata.csv'))
write.csv(tsl, paste0(procdtot, 'time_series_all.csv'))
}
#read this in
if(!overwrite)
{
  tr_wd <- read.csv(paste0(procdtot, 'fulldata.csv'))
  tsl   <- read.csv(paste0(procdtot, 'time_series_all.csv'))
}

tr_wd$ppn <- ifelse(tr_wd$ppn%in%c(41, 42), 4, tr_wd$ppn)
#tsl <- dplyr::rename(tsl, rectus_abdominis = rectus_abdominus)

Supporting results

Descriptive results

Fig. S4 shows an overview of the muscle activity patterns for each movement condition, split over weight conditions. Table S1 provides the numerical information. Table S3 provides normalized peak muscle activity by weight. All in all, these descriptive results pattern in sensible ways. Firstly, we can confirm that indeed the focal muscles powering internal (pectoralis major) or external (infraspinatus) rotation are peaking in activity during these movement conditions.

Show code

#look into NAs first
sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='expire',] #only real trials, only exhalation
numberNAs <- sum(is.na(sub$max_amp_c_around_move))
subNA <- sub %>% 
  subset(is.na(max_amp_c_around_move))
subinf <- sub %>% 
  subset(!is.finite(max_amp_c_around_move))
maxampzero <- dplyr::filter(sub, max_amp_c_around_move < 0)
minampzero <- dplyr::filter(sub, min_amp_c_around_move > 0)

There are 11 observations that are missing in the dataset. They will be excluded, as they have NA in any column on amplitude, muscle activity, or posture. Additionally, there are 2 -Inf/Inf values in the envelope, which we filter out, too.

Show code producing the plot

library(magick)

##### This the raw code for the plot
    #we have however, done after editing to prettify and mark the plot
sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='expire',] #only real trials, only exhalation

sub <- sub %>% 
   tidyr::drop_na(max_amp_c_around_move)

sub$max_pectoral <- scale(sub$max_pectoral, center = FALSE)[,1]
sub$max_infra    <- scale(sub$max_infra, center = FALSE )[,1]
sub$max_rectus   <- scale(sub$max_rectus, center = FALSE)[,1]
sub$max_erector  <- scale(sub$max_erector, center = FALSE)[,1]
sub$max_COPc     <- scale(sub$max_COPc, center = FALSE)[,1]

#reorder values
sub$movement_condition <- factor(sub$movement_condition, levels=rev(c('no movement', 'internal rotation', 'external rotation', 'flexion', 'extension')))

peaks_pectoralis_bymovement <- ggplot(sub, aes(x = movement_condition, y=max_pectoral)) +
  geom_quasirandom(aes(color = weight_condition),alpha = 0.4, dodge.width = 0.8, cex = 1) +
  geom_boxplot(aes(group=interaction(movement_condition, weight_condition)),alpha = 0, size=0.5) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Pectoralis major (internal rotator)') +
  ylab('max sEMG activity \n (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = col_pectoralis)) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm')) +
  ylim(-0.1, 7) +
  geom_hline(yintercept=median(sub$max_pectoral[sub$movement_condition == 'no movement']), linetype='dashed') +
  coord_flip()

peaks_infraspinatus_bymovement <- ggplot(sub, aes(x = movement_condition,y=max_infra)) +
  geom_quasirandom(aes(color = weight_condition), alpha = 0.4, dodge.width = 0.8, cex = 1) +
  geom_boxplot(aes(group = interaction(movement_condition, weight_condition)), alpha = 0, size=0.5) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Infraspinatus (external rotator)') +
  ylab('max sEMG activity \n (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = col_infraspinatus)) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm')) +
  ylim(-0.1, 7) +
  geom_hline(yintercept = median(sub$max_infra[sub$movement_condition == 'no movement']), linetype='dashed') +
  coord_flip()

peaks_rectus_bymovement <- ggplot(sub, aes(x = movement_condition,y=max_rectus)) +
  geom_quasirandom(aes(color = weight_condition), alpha = 0.4, dodge.width = 0.8, cex = 1) +
  geom_boxplot(aes(group=interaction(movement_condition, weight_condition)), alpha = 0,  size=0.5) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Rectus abdominis (postural)') +
  ylab('max sEMG activity \n (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = col_rectus)) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm')) +
  ylim(-0.1, 7) +
  geom_hline(yintercept=median(sub$max_rectus[sub$movement_condition == 'no movement']), linetype='dashed') +
  coord_flip()

peaks_erector_bymovement <- ggplot(sub, aes(x = movement_condition, y=max_erector)) +
  geom_quasirandom(aes(color = weight_condition), alpha = 0.4, dodge.width = 0.8, cex = 1) +
  theme_cowplot(12) +
  geom_boxplot(aes(group=interaction(movement_condition, weight_condition)), alpha = 0, size=0.5) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Erector spinae (postural)')+
  ylab('max sEMG activity \n (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = col_erector)) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm') ) +
  ylim(-0.1, 7) +
  geom_hline(yintercept=median(sub$max_erector[sub$movement_condition == 'no movement']), linetype='dashed') +
  coord_flip()

# Center of Pressure
peaks_COP_bymovement <- ggplot(sub, aes(x = movement_condition,y= max_COPc)) +
  geom_quasirandom(aes(color = weight_condition), alpha = 0.4, dodge.width = 0.8, cex = 1) +
  theme_cowplot(12) +
  geom_boxplot(aes(group=interaction(movement_condition, weight_condition)), alpha = 0, size = 0.5) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Change in center of pressure (COPc)') +
  ylab('center of pressure (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = 'grey')) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm')) +
  geom_hline(yintercept=median(sub$max_COPc[sub$movement_condition == 'no movement']),linetype='dashed') +
  coord_flip()

Figure 4: Manipulation check movement condition and peak muscle activity. The peak muscle activity, normalized for each muscle, is shown per movement and weight condition. A) highlights the fact that the pectoralis major, is - as to be expected - most active for the internal rotation movement condition. B) Also to be expected, the infraspinatus is most active during the external rotation. Note further, that both postural muscles seem more active during extension (and secondarily flexion).

Show code producing the table

library(tidyr)
library(gmodels)
library(plyr)

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='expire',] #only real trials, only exhalation

sub$max_pectoral <- scale(sub$max_pectoral, center = FALSE)[,1]
sub$max_infra    <- scale(sub$max_infra, center = FALSE )[,1]
sub$max_rectus   <- scale(sub$max_rectus, center = FALSE)[,1]
sub$max_erector  <- scale(sub$max_erector, center = FALSE)[,1]
sub$max_COPc     <- scale(sub$max_COPc, center = FALSE)[,1]

subl <- gather(sub, key = muscle, value = peak_activity, c("max_pectoral", "max_infra", "max_rectus", "max_erector"))
subl$muscle <- revalue(subl$muscle, c("max_pectoral"="pectoralis major", "max_infra"="infraspinatus", "max_rectus"="rectus abdominis", "max_erector"="erector spinae"))

# #reorder values
# subl$movement_condition <- factor(subl$movement_condition, levels=rev(c('no movement', 'internal rotation', 'external rotation', 'flexion', 'extension')))

descriptivesbymusclesmovement <- subl %>%
  group_by(muscle, movement_condition) %>%
  #I had to filter out 44 NAs
  tidyr::drop_na(peak_activity) %>% 
  dplyr::summarize(
    mean = round(ci(peak_activity)[1],2),
    sd = round(ci(peak_activity)[4],2),
    lowCI = round(ci(peak_activity)[2],2),
    hiCI = round(ci(peak_activity)[3], 2),)

#reorder to make it more comparable to previous Figure
muscle_order <- c("pectoralis major", "infraspinatus", "rectus abdominis", "erector spinae")
movement_order <- c("no movement", "internal rotation", "external rotation", "flexion", "extension")

descriptivesbymusclesmovement <- descriptivesbymusclesmovement %>%
  mutate(
    muscle = factor(muscle, levels = muscle_order),
    movement_condition = factor(movement_condition, levels = movement_order)
  )

#arrange table based on the defined order
descriptivesbymusclesmovement <- descriptivesbymusclesmovement %>%
  arrange(muscle, movement_condition)

#rename columns
descriptivesbymusclesmovement <- descriptivesbymusclesmovement %>% 
  dplyr::rename(`Muscle` = `muscle`,
                `Movement condition` = `movement_condition`,
                `Mean` = `mean`,
                `SD` = `sd`,
                `Low 95% CI` = `lowCI`,
                `High 95% CI` = `hiCI`)

Table 1: Table 2: Normalized peak muscle activity for the different movement conditions. These are the numerical results associated with Fig. S4. The values indicate the peak EMG activity normalized for each muscle.

Muscle	Movement condition	Mean	SD	Low 95% CI	High 95% CI
pectoralis major	no movement	0.19	0.02	0.16	0.22
pectoralis major	internal rotation	1.17	0.08	1.01	1.32
pectoralis major	external rotation	0.62	0.03	0.55	0.69
pectoralis major	flexion	0.85	0.05	0.75	0.95
pectoralis major	extension	0.97	0.05	0.88	1.07
infraspinatus	no movement	0.04	0.00	0.03	0.05
infraspinatus	internal rotation	0.63	0.03	0.57	0.69
infraspinatus	external rotation	1.81	0.08	1.65	1.97
infraspinatus	flexion	0.38	0.02	0.34	0.42
infraspinatus	extension	0.39	0.02	0.35	0.43
rectus abdominis	no movement	0.81	0.05	0.72	0.90
rectus abdominis	internal rotation	0.85	0.04	0.77	0.93
rectus abdominis	external rotation	0.83	0.05	0.74	0.92
rectus abdominis	flexion	0.86	0.06	0.75	0.97
rectus abdominis	extension	0.86	0.05	0.76	0.96
erector spinae	no movement	0.39	0.02	0.34	0.43
erector spinae	internal rotation	0.55	0.03	0.50	0.60
erector spinae	external rotation	0.53	0.03	0.46	0.59
erector spinae	flexion	0.90	0.07	0.77	1.04
erector spinae	extension	1.31	0.08	1.15	1.48

Further, postural muscles such as the rectus abdominis are especially active for extension movements (and secondarily flexion). Confirming their combined postural role, the muscle activity of the rectus abdominis and erector spinae are reliably correlated. Confirming their antagonistic role in rotating the humerus, the pectoralis and infraspinatus muscle activity are reliably correlated, indicative of their joint agonist/antagonistic control of posture and upper limb rotation, respectively. Lastly, from Table S3 it is clear that adding a wrist weight increases the peak muscle activity for all muscles, except the rectus abdominis (with no or barely overlapping 95 % confidence intervals in the weight vs. no weight condition).

Show code producing the table

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='expire',] #only real trials, only exhalation

sub$max_pectoral <- scale(sub$max_pectoral, center = FALSE)[,1]
sub$max_infra    <- scale(sub$max_infra, center = FALSE )[,1]
sub$max_rectus   <- scale(sub$max_rectus, center = FALSE)[,1]
sub$max_erector  <- scale(sub$max_erector, center = FALSE)[,1]
sub$max_COPc     <- scale(sub$max_COPc, center = FALSE)[,1]

subl <- gather(sub, key = muscle, value = peak_activity, c("max_pectoral", "max_infra", "max_rectus", "max_erector"))
subl$muscle <- revalue(subl$muscle, c("max_pectoral"="pectoralis major", "max_infra"="infraspinatus", "max_rectus"="rectus abdominis", "max_erector"="erector spinae"))
subl$weight_condition <- revalue(subl$weight_condition, c("no_weight"="no weight", "weight"="weight"))

#important: I had to filter out 44 NAs from peak activity or it would not work
descriptivesbymusclesweight <- subl %>%
  group_by(muscle, weight_condition) %>%
  tidyr::drop_na(peak_activity) %>% 
  dplyr::summarize(
    mean = round(ci(peak_activity)[1], 2),
    lowCI = round(ci(peak_activity)[2], 2),
    hiCI = round(ci(peak_activity)[3], 2),
    sd = round(ci(peak_activity)[4], 2))

#reorder to make it more comparable to previous Figure
muscle_order <- c("pectoralis major", "infraspinatus", "rectus abdominis", "erector spinae")
weight_order <- c("no weight", "weight")

descriptivesbymusclesweight <- descriptivesbymusclesweight %>%
  mutate(
    muscle = factor(muscle, levels = muscle_order),
    weight_condition = factor(weight_condition, levels = weight_order)
  )

#arrange table based on the defined order
descriptivesbymusclesweight <- descriptivesbymusclesweight %>%
  arrange(muscle, weight_condition)
            
#rename columns
descriptivesbymusclesweight <- descriptivesbymusclesweight %>% 
  dplyr::rename(`Muscle` = `muscle`,
                `Weight condition` = `weight_condition`,
                `Mean` = `mean`,
                `SD` = `sd`,
                `Low 95% CI` = `lowCI`,
                `High 95% CI` = `hiCI`)

#rearrange columns
descriptivesbymusclesweight <- descriptivesbymusclesweight %>%
  relocate(SD, .after = Mean)

Table 3: Table 4: Normalized peak muscle activity per weight condition. The values indicate the peak EMG activity (normalized per muscle) per weight condition.

Muscle	Weight condition	Mean	SD	Low 95% CI	High 95% CI
pectoralis major	no weight	0.64	0.03	0.58	0.70
pectoralis major	weight	0.87	0.04	0.80	0.95
infraspinatus	no weight	0.62	0.04	0.54	0.70
infraspinatus	weight	0.68	0.04	0.59	0.76
rectus abdominis	no weight	0.84	0.03	0.78	0.90
rectus abdominis	weight	0.84	0.03	0.78	0.90
erector spinae	no weight	0.60	0.03	0.53	0.66
erector spinae	weight	0.87	0.04	0.78	0.95

Exploratory visualization

Fig. S5 shows the amplitude envelope and EMG activity of a 1-second interval around the movement onset. This plot is informative about the rationale of the confirmatory analyses. Firstly, it can be seen that in general there is a positive peak during a movement onset in the amplitude envelope (relative to the trend line for that expiration), that may or may not be preceded by a less extreme negative peak. Especially, for the internal rotation we observe such a negative peak. We will therefore assess whether particular movement conditions predict higher magnitude negative and positive peaks in the amplitude envelope (analysis 1). We further assess whether particular muscles predict the magnitude of positive or negative peaks (analysis 2). Finally, we will assess whether particular muscles are related to postural stability, to confirm the posture-mediating role of the muscles (analysis 3). While we focus in this pre-registration and pilot analyses report on the peak analyses (in part to support a more straightforward power analysis), in the confirmatory study we may perform continuous trajectory analysis using generalized additive modeling similar to previous work (Pearson and Pouw 2022).

Show code producing the plot

#we will only look at exhalations and we exclude practice trials, which have a number < 9
sub2 <- subset(tsl, vocal_condition == 'expire' & trial_number > 9)
sub2$movement_condition <- revalue(sub2$movement_condition, c(
  "extension_stop" = "extension",
  "internal_rotation_stop" = "internal r.",
  "flexion_stop" = "flexion",
  "external_rotation_stop" = "external r."))

#get acceleration
sub2$unique_trial <- paste0(sub2$pp, sub2$trial_number)
sub2$speed <- scale(sub2$speed, center = FALSE)
sub2$acc <- ave(sub2$speed, sub2$unique_trial, FUN = function(x){c(0, diff(x))})


#normalize
sub2$pectoralis_major  <- scale(sub2$pectoralis_major)
sub2$infraspinatus     <- scale(sub2$infraspinatus)
sub2$rectus_abdominis  <- scale(sub2$rectus_abdominis)
sub2$erector_spinae    <- scale(sub2$erector_spinae)


#what is the average movement start and end?
average_mov_start <- mean(sub2$time_ms[sub2$movstart==0])
#this we will use to create 


smoothenvelope <- ggplot(sub2, aes(x=movstart)) +
  geom_smooth(aes(y = env_z), method = 'gam', color = 'purple') +
  facet_grid(.~movement_condition)+geom_hline(yintercept = 0) +
  xlim(-500, 500) +
  theme_cowplot(12) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('amplitude envelope') +
  xlab('time centered to movement onset (in ms)')

smoothmovementspeed <- ggplot(sub2, aes(x=movstart)) +
  geom_smooth(aes(y = speed), method = 'gam', color = 'blue') +
  facet_grid(.~movement_condition)+geom_hline(yintercept = 0) +
  xlim(-500, 500) +
  theme_cowplot(12) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('movement speed') +
  xlab('time centered to movement onset (in ms)')

smoothmovementacceleration <- ggplot(sub2, aes(x=movstart)) +
  geom_smooth(aes(y = acc), method = 'gam', color = 'red') +
  facet_grid(.~movement_condition)+geom_hline(yintercept = 0) +
  xlim(-500, 500) +
  theme_cowplot(12) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('movement acceleration') +
  xlab('time centered to movement onset (in ms)')

smoothEMGactivity <- ggplot(sub2, aes(x = movstart)) +
  geom_smooth(aes(y = pectoralis_major, color = 'pectoralis major'), method = 'gam') +
  geom_smooth(aes(y = rectus_abdominis, color = 'rectus abdominis'), method = 'gam') +
  geom_smooth(aes(y = infraspinatus, color = 'infraspinatus'), method = 'gam') +
  geom_smooth(aes(y = erector_spinae, color = 'erector spinae'), method = 'gam') +
  xlim(-500, 500) +
  facet_grid(.~movement_condition) +
  theme_cowplot(12) +
  scale_color_manual(values = colors_mus) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('EMG activity \n (normalized)') +
  xlab('time centered to movement onset (in ms)') +
  theme(legend.position = "none")    

smoothspeedCOP <- ggplot(sub2, aes(x = movstart)) +
  geom_smooth(aes(y = speed), color= 'orange', method = 'gam') +
  xlim(-500, 500) +
  facet_grid(.~movement_condition) +
  theme_cowplot(12) +
  scale_color_manual(values = colors_mus) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('speed \n (change in COP)') +
  xlab('time centered to movement onset (in ms)')

smoothCOP <- ggplot(sub2, aes(x = movstart)) +
  geom_smooth(aes(y = COPc), method = 'gam', color = 'darkgrey') +
  xlim(-500, 500) +
  facet_grid(.~movement_condition) +
  theme_cowplot(12) +
  scale_color_manual(values = colors_mus) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('change in center of pressure') +
  xlab('time centered to movement onset (in ms)')

Figure 5: Smoothed conditional means for expiration and muscle activity (with gam). Smoothed conditional means over time, where time is centered at 0 at the movement onset (as determined by the motion tracking of the wrist). The smoothed conditional means are generated by fitting non-linear smooths using generalized additive models in R-package ggplot2. It can be seen that especially for the extension movement there are clear anticipatory postural muscle activations (Aruin and Latash 1995), of the rectus abdominis before the movement onset, which is then followed by postural adjustments of the erector spinae. For the flexion condition, this activation pattern is reversed as one would expect given that the impulse vector should be directed in the opposite direction.

We here recreate the previous plot expiration amplitude (top, purple) with movement speed (middle, blue) and movement acceleration (bottom, red).

Show code producing the plot

smoothenvelope2 <- ggplot(sub2, aes(x=movstart)) +
  geom_smooth(aes(y = env_z), method = 'gam', color = '#6a3d9a') +
  facet_grid(.~movement_condition) +
  geom_hline(yintercept = 0) +
  xlim(-500, 500) +
  theme_cowplot(12) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1),
        axis.title.x = element_blank()) +
  ylab('amplitude \nenvelope')
  #xlab('time centered to movement onset (in ms)')
  # geom_vline(data = peak_values, aes(xintercept = peak), color = '#e31a1c', linetype = 'dashed', size = 1) +
  # geom_vline(data = lowpeak_values, aes(xintercept = peak), color = '#e31a1c', linetype = 'dashed', size = 1) +
  # geom_vline(data = peak_speed, aes(xintercept = peak), color = '#1f78b4', linetype = 'dashed', size = 1)

smoothmovementspeed2 <- ggplot(sub2, aes(x=movstart)) +
  geom_smooth(aes(y = speed), method = 'gam', color = '#1f78b4') +
  facet_grid(.~movement_condition) + 
  geom_hline(yintercept = 0) +
  xlim(-500, 500) +
  theme_cowplot(12) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1),
        axis.title.x = element_blank()) +
  ylab('movement \nspeed')
  #xlab('time centered to movement onset (in ms)')
  # geom_vline(data = peak_speed, aes(xintercept = peak), color = '#1f78b4', linetype = 'dashed', size = 1)

smoothmovementacceleration2 <- ggplot(sub2, aes(x=movstart)) +
  geom_smooth(aes(y = acc), method = 'gam', color = '#e31a1c') +
  facet_grid(.~movement_condition) +
  geom_hline(yintercept = 0) +
  xlim(-500, 500) +
  theme_cowplot(12) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('movement \nacceleration') +
  xlab('time centered to movement onset (in ms)')
  # geom_vline(data = peak_values, aes(xintercept = peak), color = '#e31a1c', linetype = 'dashed', size = 1) +
  # geom_vline(data = lowpeak_values, aes(xintercept = peak), color = '#e31a1c', linetype = 'dashed', size = 1)

Figure 6: Smoothed conditional means for the non-log-transformed exhalation amplitude (top, purple; z-normalized by participant), movement speed (middle, blue; z-normalized by participant), and movement acceleration (bottom, red; derivative of z-normalized speed) by movement condition across participants (via fitting non-linear smooths using generalized additive models in R-package ggplot2). Time zero is defined as movement onset (as determined by the motion tracking of the wrist.)

Main Confirmatory Analyses

As pre-registered, for the current confirmatory analyses we focus on the expiration condition, ignoring the expiration baseline conditions as they are of secondary interest at this point of our inquiry.

Do different movement conditions have different effects on expiration amplitude?

From inspecting the summarized trajectories, and the descriptive exploratory analyses above, we obtain that extension and external rotation of the arm seem to start with a negative peak around the onset of the movement, which is followed by a positive peak. Furthermore, we observe that all other movements primarily have a positive peak. A straightforward test of whether the amplitude envelope has positive or negative peaks is to assess differences in peaks in the expiration conditions per movement (and weight) condition.

Show code producing the plot

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='expire',] #only real trials, only exhalation

#reorder values & turn into factor
sub$movement_condition <- factor(sub$movement_condition, levels=c(
  "no movement",
  "internal rotation",
  "external rotation",
  "flexion",
  "extension"))
#rename no_weight in no weight
sub$weight_condition <- dplyr::recode(sub$weight_condition,
                               weight = "weight",
                               no_weight = "no weight")
#turn into factor
sub$weight_condition <- factor(sub$weight_condition, levels=c(
  "no weight",
  "weight"))

#remove both -Inf (=2) and NA (=11)
sub$max_amp_c_around_move <- na_if(sub$max_amp_c_around_move, -Inf)
sub <- sub %>% 
  tidyr::drop_na(max_amp_c_around_move)

#this also removes the NA showing up in min_amp so no need to filter them out, too
# sub$min_amp_c_around_move <- na_if(sub$min_amp_c_around_move, -Inf)
# sub <- sub %>% 
#   tidyr::drop_na(max_amp_c_around_move)


pospeakamplitudemovement <- ggplot(sub, aes(x = movement_condition, y = max_amp_c_around_move)) +
  geom_quasirandom(size = 2,alpha = 0.5, dodge.width=0.8, cex=1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("movement condition") +
  ggtitle("Positive peak amplitude") +
  ylab("amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title  =element_blank()) +
  geom_hline(yintercept = median(sub$max_amp_c_around_move[sub$movement_condition == "no movement"]), size=2, alpha=0.5, color = "blue")
 # ylim(-0.01, 0.35)

negpeakamplitudemovement <- ggplot(sub, aes(x = movement_condition, y = min_amp_c_around_move)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("movement condition") +
  ggtitle('Negative peak amplitude') +
  ylab("amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title = element_blank()) +
  geom_hline(yintercept = median(sub$min_amp_c_around_move[sub$movement_condition == "no movement"]), size = 2, alpha = 0.5, color = "blue")
  #ylim(-0.15, 0.01)

pospeakamplitudeweight <- ggplot(sub[sub$movement_condition != 'no movement', ], aes(x = weight_condition, y = max_amp_c_around_move)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("weight condition") +
  ggtitle("Positive peak amplitude") +
  ylab("amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title=element_blank())
  #ylim(-0.01, 0.35)

negpeakamplitudeweight <- ggplot(sub[sub$movement_condition != 'no movement', ], aes(x = weight_condition,y = min_amp_c_around_move)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("weight condition") +
  ggtitle("Negative peak amplitude") +
  ylab("amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title=element_blank())
  #ylim(-0.15, 0.01)

maxampzero <- dplyr::filter(sub, max_amp_c_around_move < 0)
minampzero <- dplyr::filter(sub, min_amp_c_around_move > 0)

Figure 7: Effects of movement condition of positive and negative peaks in the voice amplitude. The upper part shows the positive peaks in the amplitude envelope during the different movement and weight conditions. The lower part shows the negative peaks (hence the negative values; note, in the modeling we will absolutize these values). It is clear that relative to expiration of no movement, there are especially positive, but also more negative peaks in the amplitude envelope for the different movement conditions. Please not that the y-axis is different for negative peaks because the magnitude of them is lower than for positive peaks.

Since most of the data points are very close to 0 and the distributions have long tails, we will employ a log-transformation. For this, we first filter out the 50 values in the positive peaks that are below zero. We also filter out the 3 values in the negative peaks that are above zero, before we absolutize the numbers. Then, the log-transformation is employed. This is done, as we cannot log-transform negative values and simply absolutizing them all would inflate potential effects.

Show code for log-transform and producing the plot

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='expire',] #only real trials, only exhalation

#reorder values & turn into factor
sub$movement_condition <- factor(sub$movement_condition, levels=c(
  "no movement",
  "internal rotation",
  "external rotation",
  "flexion",
  "extension"))

#rename them for plot
# sub$movement_condition <- dplyr::recode(sub$movement_condition,
#                                `no movement` = "no move.",
#                                `internal rotation` = "int. rot",
#                                `external rotation` = "ext. rot" ,
#                                `flexion` = "flexion",
#                                `extension` = "extension")

#rename no_weight in no weight
sub$weight_condition <- dplyr::recode(sub$weight_condition,
                               weight = "weight",
                               no_weight = "no weight")
#turn into factor
sub$weight_condition <- factor(sub$weight_condition, levels=c(
  "no weight",
  "weight"))

#remove both -Inf (=2) and NA (=11)
sub$max_amp_c_around_move <- na_if(sub$max_amp_c_around_move, -Inf)
sub <- sub %>% 
  tidyr::drop_na(max_amp_c_around_move)

#filter out values in max_amp_c_around_move that are below 0 and then log-transform
submax <- sub %>% 
  dplyr::filter(max_amp_c_around_move > 0)
submax$max_amp_log <- log(submax$max_amp_c_around_move)

#filter out values in min_amp_c_around_move that are above 0, absolutize, and then log-transform
submin <- sub %>% 
  dplyr::filter(min_amp_c_around_move < 0)
submin$min_amp_log <- log(abs(submin$min_amp_c_around_move))

logpospeakamplitudemovement <- ggplot(submax, aes(x = movement_condition, y = max_amp_log)) +
  geom_quasirandom(size = 2,alpha = 0.5, dodge.width=0.8, cex=1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("movement condition") +
  ggtitle("Positive peak amplitude") +
  ylab("log-transformed amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title =element_blank()) +
  geom_hline(yintercept = median(submax$max_amp_log[submax$movement_condition == "no move."]), size=2, alpha=0.5, color = "blue")
 # ylim(-0.01, 0.35)

lognegpeakamplitudemovement <- ggplot(submin, aes(x = movement_condition, y = min_amp_log)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("movement condition") +
  ggtitle('Negative peak amplitude') +
  ylab("log-transformed amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title = element_blank()) +
  geom_hline(yintercept = median(submin$min_amp_log[submin$movement_condition == "no move."]), size = 2, alpha = 0.5, color = "blue") 
  #ylim(-0.15, 0.01)

logpospeakamplitudeweight <- ggplot(submax[submax$movement_condition != 'no move.', ], aes(x = weight_condition, y = max_amp_log)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("weight condition") +
  ggtitle("Positive peak amplitude") +
  ylab("log-transformed amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title=element_blank())
  #ylim(-0.01, 0.35)

lognegpeakamplitudeweight <- ggplot(submin[submin$movement_condition != 'no move.', ], aes(x = weight_condition,y = min_amp_log)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("weight condition") +
  ggtitle("Negative peak amplitude") +
  ylab("log-transformed amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title=element_blank()) 
  #ylim(-0.15, 0.01)

Figure 8: Effects of movement condition of positive and negative peaks in the voice amplitude. Same as the previous plot, but we here log-transformed the amplitude peak values.

Show code running models on log-transformed positive amplitude peaks

#sub <- tr_wd[tr_wd$trialindex>9 & tr_wd$vocal_condition=='expire', ] #only real trials, only exhalation
#sub$movement_condition <- factor(sub$movement_condition, levels=c('no movement', 'internal rotation', 'external rotation', 'flexion', 'extension'))
submax$weight_condition <- dplyr::recode(submax$weight_condition,
                               weight = "weight",
                               no_weight = "no weight")
submax$weight_condition <- factor(submax$weight_condition, levels=c("no weight", "weight"))

submin$weight_condition <- dplyr::recode(submin$weight_condition,
                               weight = "weight",
                               no_weight = "no weight")
submin$weight_condition <- factor(submin$weight_condition, levels=c("no weight", "weight"))

#the Inf values should already be filtered out, this acts as a safety net
# sub$max_amp_log <- na_if(sub$max_amp_log, Inf)
# sub <- sub %>%
#   tidyr::drop_na(max_amp_log)


#model 1

pospeaks_nullmodel <- lmer(max_amp_log ~ 1 + (1|ppn), data = submax)
pospeaks_model   <- lmer(max_amp_log ~ weight_condition + movement_condition + (1|ppn), data = submax)
pospeaks_interactionmodel   <- lmer(max_amp_log ~ weight_condition * movement_condition + (1|ppn), data = submax)
# pospeaks_nullmodel <- lmer(max_amp_c_around_move ~ 1+(1|ppn), data = sub)
# pospeaks_model   <- lmer(max_amp_c_around_move ~ weight_condition * movement_condition + (1|ppn), data = sub)

# in the vocalize condition, this used to be an additive model of both movement & weight
#logpospeaks_model   <- lmer(max_amp_log ~ weight_condition + movement_condition + (1|ppn), data = sub)
#running "weight_condition * movement_condition" returns no significant interactions with weight
pospeaks_comparison <- anova(pospeaks_nullmodel, pospeaks_model)


modelsummarypospeaksweightmovement <- summary(pospeaks_model)$coefficients
colnames(modelsummarypospeaksweightmovement) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarypospeaksweightmovement) <- c(
  'Intercept (no movement/no weight)',
  'vs. weight',
  'vs. internal rotation',
  'vs. external rotation',
  'vs. flexion',
  'vs. extension')

#posthoc
pospeaks_model_posthoc <- emmeans(pospeaks_model, list(pairwise ~ movement_condition), adjust = "tukey")
#pcomparisons <- posthoc$`pairwise differences of movement_condition`

  #########################################also make a comparison with only movement conditions for the weight condition (CHANGECONF)
subsubmax <- subset(submax, movement_condition != 'no movement')
#model 1
pospeaks_nullmodel_onlymovement <- lmer(max_amp_log ~ 1 + (1 | ppn), data = subsubmax)
pospeaks_model_onlymovement   <- lmer(max_amp_log ~ weight_condition + (1 | ppn), data = subsubmax)
#now, it's actually weight vs no weight that produces a better model than null
#but none of the factors are significant, either in movement model, weight model, or interaction model
pospeaks_comparison_onlymovement <- anova(pospeaks_nullmodel_onlymovement, pospeaks_model_onlymovement)

#table matrix
  #Estimate, 

modelsummarypospeaksweightmovement_onlymovement <- summary(pospeaks_model_onlymovement)$coefficients
colnames(modelsummarypospeaksweightmovement_onlymovement) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarypospeaksweightmovement_onlymovement) <- c('Intercept (no weight)', 'vs. weight')
#`r printnum(ifelse(pospeaks_comparison[2,8] > .001, paste('=', pospeaks_comparison[2,8]), '< .001'))`)

We first modeled with a mixed linear regression the variation in positive peaks in the amplitude envelope with the predictors movement condition and weight condition (using R-package lme4), with participant as random intercept (for more complex random slope models did not converge). A model with the movement conditions explained more variance than a base model predicting the overall mean, Change in Chi-Squared (5.00) = 90.67, p = < .001. Trying interactions in the model between weight and movement conditions did not lead to a significant main effect of weight or significant interactions between the two factors. The model coefficients are given in Table S5. There is thus no effect of of wrist weight in this sample. Further, all movements (extension, flexion, internal rotation, external rotation) lead to statistically reliable increases in positive peaks in the amplitude envelope relative to the no movement condition.

Table 5: Table 6: Effects of weight and movement condition on magnitude positive peaks in amplitude envelope.

	Estimate	SE	df	t-value	p-value
Intercept (no movement/no weight)	-5.791	0.329	20.339	-17.600	0.000
vs. weight	0.114	0.099	553.313	1.153	0.249
vs. internal rotation	1.313	0.156	552.895	8.395	0.000
vs. external rotation	1.096	0.158	553.076	6.933	0.000
vs. flexion	1.283	0.156	552.919	8.202	0.000
vs. extension	1.113	0.156	552.994	7.152	0.000

Note that in the previous model weight did not turn out a significant factor. During the no movement condition the weight is not affecting the participant. Thus, a better estimate of the effect of wrist weight is to assess weight vs. no weight conditions for only the movement conditions, thus excluding the no movemement condition. However, Table S7 below shows that this comparison does not lead to weight being a significant factor either.

Table 7: Table 8: Effects of weight on magnitude of positive peaks in amplitude envelope (no movement condition filtered out).

	Estimate	SE	df	t-value	p-value
Intercept (no weight)	-4.579	0.293	16.815	-15.633	0.000
vs. weight	0.129	0.099	445.220	1.301	0.194

Show code running models on log-transformed negative amplitude peaks

# model 2
# sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='expire',] #only real trials, only exhalation
# sub$movement_condition <- factor(sub$movement_condition, levels=c('no movement', 'internal rotation', 'external rotation', 'flexion', 'extension'))
submin$weight_condition <- dplyr::recode(submin$weight_condition,
                               weight = "weight",
                               no_weight = "no weight")
submin$weight_condition <- factor(submin$weight_condition, levels=c("no weight", "weight"))

#sub$min_amp_c_around_move_abs <- abs(sub$min_amp_c_around_move)

negpeaks_nullmodel <- lmer(min_amp_log ~  1  + (1|ppn), data = submin)
negpeaks_model <- lmer(min_amp_log ~ weight_condition + movement_condition + (1|ppn), data = submin)
#negpeaks_interactionmodel   <- lmer(min_amp_log ~ weight_condition * movement_condition + (1|ppn), data = submin)
#again, it's only movement that has significant factors (vs weight or significant interactions in weight * movement)
negpeaks_comparison <- anova(negpeaks_nullmodel, negpeaks_model)


#table matrix
modelsummarynegpeaksweightmovement <- summary(negpeaks_model)$coefficients
colnames(modelsummarynegpeaksweightmovement) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarynegpeaksweightmovement) <- c(
'Intercept (no movement/no weight)',
'vs. weight',
'vs. internal rotation',
'vs. external rotation',
'vs. flexion',
'vs. extension')

#also make a comparison with only movement conditions for the weight condition (CHANGECONF)
subsubmin <- subset(submin, movement_condition != 'no movement')

negpeaks_nullmodel_onlymovement <- lmer(min_amp_log ~ 1 + (1 | ppn), data = subsubmin)
negpeaks_model_onlymovement   <- lmer(min_amp_log ~ weight_condition + (1 | ppn), data = subsubmin)
negpeaks_comparison_onlymovement <- anova(negpeaks_nullmodel_onlymovement, negpeaks_model_onlymovement)

modelsummarynegpeaksweightmovement_onlymovement <- summary(negpeaks_model_onlymovement)$coefficients
colnames(modelsummarynegpeaksweightmovement_onlymovement) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarynegpeaksweightmovement_onlymovement) <- c('Intercept (no weight)', 'vs. weight')

Secondly, we modeled in a similar way the negative peaks in the amplitude envelope, and found that neither a model with movement conditions, weight and movement conditions, or interactions was able to explain more variance than a base model predicting the overall mean, for the additive model Change in Chi-Squared (5.00) = 71.83, p = < .001. Still, the model coefficients are shown in Table S9.

Table 9: Table 10: Effects of weight and movement condition on magnitude negative peaks in amplitude envelope.

	Estimate	SE	df	t-value	p-value
Intercept (no movement/no weight)	-6.526	0.443	17.828	-14.743	0.000
vs. weight	0.123	0.091	600.217	1.354	0.176
vs. internal rotation	0.926	0.143	599.973	6.479	0.000
vs. external rotation	0.699	0.144	599.976	4.866	0.000
vs. flexion	1.038	0.143	599.977	7.259	0.000
vs. extension	1.036	0.144	599.978	7.199	0.000

As above, we tried isolating weight (vs no weight) as a factor in only movement conditions (i.e. no movement filtered out). However, as above for positive amplitude peaks, we did not find a significant effect of weight on the magnitude of negative peaks.

Table 11: Table 12: Effects of weight on magnitude of negative peaks in amplitude envelope (no movement condition filtered out).

	Estimate	SE	df	t-value	p-value
Intercept (no weight)	-5.580	0.403	16.406	-13.832	0.00
vs. weight	0.108	0.090	481.163	1.202	0.23

Is different muscle activity differently affecting amplitude of expiration?

Show code calculating VIF

sub <- tr_wd[tr_wd$vocal_condition == 'expire' & tr_wd$trialindex>9,] #only real trials, only exhalation
#check VIF
library(usdm)
df = cbind.data.frame(sub$max_COPc, sub$max_erector, sub$max_infra, sub$max_pectoral, sub$max_rectus) # Data Frame
VIFFIES <- vif(df)

Since each movement and weight condition is designed to recruit different muscle activations, we can also directly relate muscle activity peaks with the positive and negative peaks in the amplitude envelope. We use a similar linear mixed regression approach to model variance in vocal amplitude peaks with peaks in muscle activity for the different muscles measured. Since there are correlations between these muscle activities (see Tables S19, S21, S23, S25), we first assessed the VIF’s between the muscle activity peaks, which yielded a maximum VIF value of 1.25. Since this is considered a low value (VIF > 5 is generally considered problematic), we can combine the different muscle activity measurements in one model to predict amplitude envelope peaks. Figure S9 shows the graphical results of these relationships.

Show code producing the plot

#we here recycle submax & submin

#normalize submax
submax$max_pectoral  <- scale(submax$max_pectoral, center = FALSE)
submax$max_infra     <- scale(submax$max_infra, center = FALSE)
submax$max_rectus    <- scale(submax$max_rectus, center = FALSE)
submax$max_erector   <- scale(submax$max_erector, center = FALSE)

#normalize submin
submin$max_pectoral  <- scale(submin$max_pectoral, center = FALSE)
submin$max_infra     <- scale(submin$max_infra, center = FALSE)
submin$max_rectus    <- scale(submin$max_rectus, center = FALSE)
submin$max_erector   <- scale(submin$max_erector, center = FALSE)


#remove both -Inf (=2) and NA (=11)
# submax$maxpectoral <- na_if(submax$maxpectoral, -Inf)
# submax <- submax %>% 
#   tidyr::drop_na(maxpectoral)

#filter out values in max_amp_c_around_move that are below 0 and then log-transform
submaxpositive <- submax %>% 
  dplyr::filter(max_pectoral > 0)
submaxpositive <- submaxpositive %>% 
  dplyr::filter(max_infra > 0)
submaxpositive <- submaxpositive %>% 
  dplyr::filter(max_erector > 0)
submaxpositive <- submaxpositive %>% 
  dplyr::filter(max_rectus > 0)

submaxpositive$max_pectoral_log <- log(submaxpositive$max_pectoral)
submaxpositive$max_infra_log <- log(submaxpositive$max_infra)
submaxpositive$max_erector_log <- log(submaxpositive$max_erector)
submaxpositive$max_rectus_log <- log(submaxpositive$max_rectus)

#filter out values in min_amp_c_around_move that are above 0, absolutize, and then log-transform
submin <- sub %>% 
  dplyr::filter(min_amp_c_around_move < 0)
submin$min_amp_log <- log(abs(submin$min_amp_c_around_move))



pospeakamplitude_peakpectoralis <- ggplot(submax, aes(y = max_amp_log)) +
  geom_point(aes(x = max_pectoral, color = 'pectoralis major'), size = 0.2) +
  geom_smooth(aes(x = max_pectoral, color = 'pectoralis major'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('positive peak \n amplitude envelope') +
  xlab('peak sEMG \n pectoralis major') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7) +
  ylim(-8.5, 1.5)

pospeakamplitude_peakrectus <- ggplot(submax, aes(y = max_amp_log)) +
  geom_point(aes(x = max_rectus, color = 'rectus abdominis'), size = 0.2) +
  geom_smooth(aes(x = max_rectus, color = 'rectus abdominis'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('positive peak \n amplitude envelope') +
  xlab('peak sEMG \n rectus abdominis') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7) +
  ylim(-8.5, 1.5)

pospeakamplitude_peakerector <- ggplot(submax, aes(y = max_amp_log)) +
  geom_point(aes(x = max_erector, color = 'erector spinae'), size = 0.2) +
  geom_smooth(aes(x = max_erector, color = 'erector spinae'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('positive peak \n amplitude envelope') +
  xlab('peak sEMG \n erector spinae') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7) +
  ylim(-8.5, 1.5)

pospeakamplitude_peakinfraspinatus <- ggplot(submax, aes(y = max_amp_log)) +
  geom_point(aes(x = max_infra, color = 'infraspinatus'), size = 0.2) +
  geom_smooth(aes(x = max_infra, color = 'infraspinatus'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('positive peak \n amplitude envelope') +
  xlab('peak sEMG \n infraspinatus') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7) +
  ylim(-8.5, 1.5)


negpeakamplitude_peakpectoralis <- ggplot(submin, aes(y = min_amp_log)) +
  geom_point(aes(x = max_pectoral, color = 'pectoralis major'), size = 0.2) +
  geom_smooth(aes(x = max_pectoral, color = 'pectoralis major'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('negative peak \n amplitude envelope') +
  xlab('peak sEMG \n pectoralis major') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7) +
  ylim(-8.5, 1.5)

negpeakamplitude_peakrectus <- ggplot(submin, aes(y = min_amp_log)) +
  geom_point(aes(x = max_rectus, color = 'rectus abdominis'), size = 0.2) +
  geom_smooth(aes(x = max_rectus, color = 'rectus abdominis'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position="none") +
  ylab('negative peak \n amplitude envelope') +
  xlab('peak sEMG \n rectus abdominis') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7) +
  ylim(-8.5, 1.5)

negpeakamplitude_peakerector <- ggplot(submin, aes(y = min_amp_log)) +
  geom_point(aes(x = max_erector, color = 'erector spinae'), size = 0.2) +
  geom_smooth(aes(x = max_erector, color = 'erector spinae'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('negative peak \n amplitude envelope') +
  xlab('peak sEMG \n erector spinae') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7) +
  ylim(-8.5, 1.5)

negpeakamplitude_peakeinfraspinatus <-ggplot(submin, aes(y = min_amp_log)) +
  geom_point(aes(x = max_infra, color = 'infraspinatus'), size =0.2) +
  geom_smooth(aes(x = max_infra, color = 'infraspinatus'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('negative peak \n amplitude envelope') +
  xlab('peak sEMG \n infraspinatus') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7) +
  ylim(-8.5, 1.5)

Figure 9: Relations between peak muscle activity and positive peaks in the amplitude envelope

Show code modeling positive amplitude peaks by muscle

# sub <- tr_wd[tr_wd$vocal_condition == 'expire' & tr_wd$trialindex>9, ] #only real trials, only exhalation
# 
# #onyl needed bc max_amp_c_around_move uses envdb
# # sub$max_amp_c_around_move <- na_if(sub$max_amp_c_around_move, -Inf)
# sub <- sub %>% 
#   tidyr::drop_na(max_amp_c_around_move)

#recycle submax & submin again

#model 1
pospeaks_bymuscle_nullmodel <- lmer(max_amp_log ~ 1 + (1|ppn), data = submax)
pospeaks_bymuscle <- lmer(max_amp_log ~ max_erector + max_infra + max_pectoral + max_rectus + (1 | ppn), data = submax)
pospeaks_bymuscle_comparison <- anova(pospeaks_bymuscle_nullmodel, pospeaks_bymuscle)

#get model summary
modelsummarypospeaksbymuscle <- summary(pospeaks_bymuscle)$coefficients
colnames(modelsummarypospeaksbymuscle) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarypospeaksbymuscle) <- c('Intercept', 'erector spinae', 'infraspinatus', 'pectoralis major', 'rectus abdominis')

A model with participant as random intercept (a more complex random slope model did not converge) the different peak muscle activities explained more variance than a base model predicting the overall mean, Change in Chi-Squared (4.00) = 51.87, p = < .001). The model coefficients are given in Table S13. Further, Table S13 shows that peak EMG activity in the infraspinatus leads to statistically reliable increases in positive peaks in the amplitude envelope (with stronger effects).

Table 13: Table 14: Linear mixed regression model assessing the relation of peak muscle activity with the positive peak in the amplitude envelope.

	Estimate	SE	df	t-value	p-value
Intercept	-5.430	0.337	21.580	-16.124	0.000
erector spinae	0.186	0.083	557.812	2.232	0.026
infraspinatus	0.276	0.069	555.346	3.995	0.000
pectoralis major	0.421	0.083	555.750	5.047	0.000
rectus abdominis	0.030	0.113	563.642	0.264	0.792

We similarly modeled the variance of the magnitude in negative peaks in the amplitude envelope with muscle activity peaks (see right side of Figure S9 for graphical results), and observed that such a model did not perform better than a base model predicting the overall mean, Change in Chi-Squared = 40.40 (4.00), p = < .001. The model coefficients are given in Table S15.

Table 15: Table 16: Linear mixed regression model assessing the relation of peak muscle activity with the negative peak in the amplitude envelope. The negative peaks in the amplitude envelope are absolutized.

	Estimate	SE	df	t-value	p-value
Intercept	-6.275	0.452	18.398	-13.881	0.000
erector spinae	0.043	0.019	603.162	2.307	0.021
infraspinatus	0.018	0.011	601.662	1.618	0.106
pectoralis major	0.070	0.014	601.935	5.068	0.000
rectus abdominis	0.019	0.030	605.736	0.633	0.527

Is a particular muscle related to postural stability?

Show code modeling change in posture by muscle

sub <- tr_wd[tr_wd$vocal_condition == 'expire' & tr_wd$trialindex>9, ] #only real trials, only exhalation
sub$max_COPc <- scale(sub$max_COPc)
sub$max_pectoral  <- scale(sub$max_pectoral, center = FALSE)
sub$max_infra     <- scale(sub$max_infra, center = FALSE)
sub$max_rectus    <- scale(sub$max_rectus, center = FALSE)
sub$max_erector   <- scale(sub$max_erector, center = FALSE)

maxCOPc_bymuscle_nullmodel <- lmer(max_COPc ~ 1 + (1 | ppn), data = sub)
maxCOPc_bymuscle <- lmer(max_COPc ~ max_erector + max_infra + max_pectoral + max_rectus + (1 | ppn), data = sub)
maxCOPc_bymuscle_comparison <- anova(maxCOPc_bymuscle_nullmodel, maxCOPc_bymuscle)

#get model summary and prepare table
modelsummarymaxCOPcbymuscle <- summary(maxCOPc_bymuscle)$coefficients
colnames(modelsummarymaxCOPcbymuscle) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarymaxCOPcbymuscle) <- c('Intercept', 'erector spinae', 'infraspinatus', 'pectoralis major', 'rectus abdominis')

Finally, we will assess which muscle activity can be related to changes in the center of pressure, which would directly confirm that gesture-speech biomechanics relate to postural stability (Pouw, Harrison, and Dixon 2020). Figure S10 shows the graphical results. We similarly performed a linear mixed regression (with participant as random intercept) with a model containing peak EMG activity for each muscle which was regressed on the peak in change in the center of pressure (COPc). We obtained that a base model predicting the over all mean of COPc was outperformed relative to said model, Change in Chi-Squared (4.00) = 153.38, p < .001). Table S17 provides the coefficient information. We find that the erector spinae, infraspinatus, and pectoralis major indeed reliably predict the magnitude of changes in the center of pressure, while the rectus abdominis does not reliably relate to the changes in ground reaction forces.

Show code producing the plot

sub <- tr_wd[tr_wd$trialindex>9,] #only real trials
sub <- subset(sub, movement_condition != 'no movement' & vocal_condition == 'expire') #only movement, only exhalation

sub <- sub %>% 
   tidyr::drop_na(max_amp_c_around_move)

sub$max_pectoral <- scale(sub$max_pectoral)
sub$max_infra    <- scale(sub$max_infra )
sub$max_rectus   <- scale(sub$max_rectus)
sub$max_erector  <- scale(sub$max_erector)

#correlation plot, to confirm role of postural muscle activations
COPcbymuscle_infraspinatus <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_infra, color = 'infraspinatus'), size = 0.2) +
  geom_smooth(aes(y = max_infra, color = 'infraspinatus'), method='lm') +
  ggtitle('Focal muscles: \n infraspinatus') +
  scale_color_manual(values = colors_mus) +
  ylab('max sEMG activity \n (normalized)') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 12, color = col_infraspinatus),
        axis.text.x = element_text(size = 8),
        legend.position = "none")

COPcbymuscle_pectoralis <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_pectoral, color = 'pectoralis major'), size=0.2) +
  geom_smooth(aes(y = max_pectoral, color = 'pectoralis major'), method = 'lm') +
  ggtitle('Focal muscles: \n pectoralis major') +
  scale_color_manual(values = colors_mus) +
  ylab('max sEMG activity \n (normalized)') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 12, color = col_pectoralis),
        axis.text.x = element_text(size = 8),
        legend.position = "none")

COPcbymuscle_rectus <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_rectus, color = 'rectus abdominis'), size = 0.2) +
  geom_smooth(aes(y = max_rectus, color = 'rectus abdominis'), method = 'lm') +
  #ggtitle('Postural muscles') +
  scale_color_manual(values = colors_mus) +
  ylab('max sEMG activity \n (normalized)') +
  ggtitle('Postural muscles: \n rectus abdominis') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 12, color = col_rectus),
        axis.text.x = element_text(size = 8),
        legend.position = "none")


COPcbymuscle_erector <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_erector, color = 'erector spinae'), size = 0.2) +
  geom_smooth(aes(y = max_erector, color = 'erector spinae'), method = 'lm') +
  #ggtitle('Postural muscles') +
  scale_color_manual(values = colors_mus) +
  ylab('max sEMG activity \n (normalized)') +
  ggtitle('Postural muscles: \n erector spinae') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 12, color = col_erector),
        axis.text.x = element_text(size = 8),
        legend.position = "none")

Figure 10: Confirmation of postural muscle activation during changes in center of mass.

Table 17: Table 18: Linear mixed regression model for predicting peak change in center of pressure based on muscle activity.

	Estimate	SE	df	t-value	p-value
Intercept	-0.748	0.122	51.152	-6.156	0.000
erector spinae	0.452	0.055	621.673	8.288	0.000
infraspinatus	0.216	0.045	613.744	4.775	0.000
pectoralis major	0.387	0.055	615.558	7.047	0.000
rectus abdominis	-0.048	0.072	604.185	-0.660	0.509

Other exploratory analyses

Along with the confirmatory analyses reported above that relate to questions posed in the pre-registration, we report on additional, exploratory analyses. Firstly we assess for possible context-dependent interrelationships of posture, muscle activity, and voice, by producing correlation matrices per movement condition (see Tables S19, S21, S23, S25). We see that producing different movements can recruit variable interrelationships between muscle, posture and voice. For exhalation, we use the log-transformed amplitude.

Show code producing the correlation matrices and plots

library(corx)


#onyl needed bc max_amp_c_around_move uses envdb
# sub$max_amp_c_around_move <- na_if(sub$max_amp_c_around_move, -Inf)
sub <- sub %>% 
  tidyr::drop_na(max_amp_c_around_move)

#correlation between positive and negative peaks in amplitude and change in center of pressure
correlation_pospeakCOPc_infraspinatus <- ggplot(submax, aes(x = max_COPc)) +
  geom_point(aes(y = max_amp_log, color = 'infraspinatus'), size = 0.2) +
  geom_smooth(aes(y = max_amp_log), method = 'lm') +
  ggtitle('Focal muscles: \n infraspinatus') +
  scale_color_manual(values = colors_mus) +
  ylab('positive amplitude peaks') +
  xlab('max. change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 20, color = col_infraspinatus),
        legend.position = "none") #+facet_grid(.~movement_condition)

correlation_negpeakCOPc_infraspinatus <- ggplot(submin, aes(x = max_COPc)) +
  geom_point(aes(y = min_amp_log, color = 'infraspinatus'), size = 0.2) +
  geom_smooth(aes(y = min_amp_log), method = 'lm') +
  ggtitle('Focal muscles: \n infraspinatus') +
  scale_color_manual(values = colors_mus) +
  ylab('negative amplitude peaks') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 20, color = col_infraspinatus),
        legend.position = "none") #+facet_grid(.~movement_condition)

#internal rotation
condition_internalrotation <- submax$movement_condition=='internal rotation'
data_condition_internalrotation <- cbind.data.frame(submax$max_pectoral[condition_internalrotation], submax$max_infra[condition_internalrotation], submax$max_rectus[condition_internalrotation], submax$max_erector[condition_internalrotation], submax$max_COPc[condition_internalrotation], submax$max_amp_log[condition_internalrotation])

#rename columns
data_condition_internalrotation_rn <- data_condition_internalrotation %>% 
  dplyr::rename(
    `max. EMG pectoralis major` = `submax$max_pectoral[condition_internalrotation]`,
    `max. EMG infraspinatus` = `submax$max_infra[condition_internalrotation]`,
    `max. EMG rectus abdominis` = `submax$max_rectus[condition_internalrotation]`,
    `max. EMG erector spinae` = `submax$max_erector[condition_internalrotation]`,
    `max. change in center of pressure` = `submax$max_COPc[condition_internalrotation]`,
    `max. positive amplitude expiration` = `submax$max_amp_log[condition_internalrotation]`
    )
    
corr_internalrotation <- data_condition_internalrotation_rn %>%
  corx::corx(triangle = "lower",
       stars = c(0.05, 0.01, 0.001),
       note = "Internal rotation") 
#knitr::kable(corr_internalrotation$apa)


#external rotation
condition_externalrotation <- submax$movement_condition=='external rotation'
data_condition_externalrotation <- cbind.data.frame(submax$max_pectoral[condition_externalrotation], submax$max_infra[condition_externalrotation], submax$max_rectus[condition_externalrotation], submax$max_erector[condition_externalrotation], submax$max_COPc[condition_externalrotation], submax$max_amp_log[condition_externalrotation])

#rename columns
data_condition_externalrotation_rn <- data_condition_externalrotation %>% 
  dplyr::rename(
    `max. EMG pectoralis major` = `submax$max_pectoral[condition_externalrotation]`,
    `max. EMG infraspinatus` = `submax$max_infra[condition_externalrotation]`,
    `max. EMG rectus abdominis` = `submax$max_rectus[condition_externalrotation]`,
    `max. EMG erector spinae` = `submax$max_erector[condition_externalrotation]`,
    `max. change in center of pressure` = `submax$max_COPc[condition_externalrotation]`,
    `max. positive amplitude expiration` = `submax$max_amp_log[condition_externalrotation]`
    )

corr_externalrotation <- data_condition_externalrotation_rn %>%
  corx::corx(triangle = "lower",
       stars    = c(0.05, 0.01, 0.001),
       note = "External rotation") 
#knitr::kable(corr_externalrotation$apa)


#extension
condition_extension <- submax$movement_condition=='extension'
data_condition_extension <- cbind.data.frame(submax$max_pectoral[condition_extension], submax$max_infra[condition_extension], submax$max_rectus[condition_extension], submax$max_erector[condition_extension], submax$max_COPc[condition_extension], submax$max_amp_log[condition_extension])
#cor.test(sub$max_rectus[condition_extension], sub$max_COPc[condition_extension])

#rename columns
data_condition_extension_rn <- data_condition_extension %>% 
  dplyr::rename(
    `max. EMG pectoralis major` = `submax$max_pectoral[condition_extension]`,
    `max. EMG infraspinatus` = `submax$max_infra[condition_extension]`,
    `max. EMG rectus abdominis` = `submax$max_rectus[condition_extension]`,
    `max. EMG erector spinae` = `submax$max_erector[condition_extension]`,
    `max. change in center of pressure` = `submax$max_COPc[condition_extension]`,
    `max. positive amplitude expiration` = `submax$max_amp_log[condition_extension]`
    )

corr_extension <- data_condition_extension_rn %>%
  corx::corx(triangle = "lower",
       stars    = c(0.05, 0.01, 0.001),
       note = "Extension")
#knitr::kable(corr_extension$apa)


#flexion
condition_flexion <- submax$movement_condition=='flexion'
data_condition_flexion <- cbind.data.frame(submax$max_pectoral[condition_flexion], submax$max_infra[condition_flexion], submax$max_rectus[condition_flexion], submax$max_erector[condition_flexion], submax$max_COPc[condition_flexion], submax$max_amp_log[condition_flexion])

#rename columns
data_condition_flexion_rn <- data_condition_flexion %>% 
  dplyr::rename(
    `max. EMG pectoralis major` = `submax$max_pectoral[condition_flexion]`,
    `max. EMG infraspinatus` = `submax$max_infra[condition_flexion]`,
    `max. EMG rectus abdominis` = `submax$max_rectus[condition_flexion]`,
    `max. EMG erector spinae` = `submax$max_erector[condition_flexion]`,
    `max. change in center of pressure` = `submax$max_COPc[condition_flexion]`,
    `max. positive amplitude expiration` = `submax$max_amp_log[condition_flexion]`
    )

corr_flexion <- data_condition_flexion_rn %>%
  corx::corx(triangle = "lower",
       stars    = c(0.05, 0.01, 0.001),
       note = "Flexion") 
#knitr::kable(corr_flexion$apa)


#very stable correlation of posture and internal rotation
correlation_posture_internalrotation <- submax %>% 
  dplyr::filter(movement_condition == "internal rotation") %>% 
  ggplot(aes(x = max_COPc, y = max_amp_log)) +
  geom_point(size = 3) +
  geom_smooth(method = 'lm', color = 'red', size = 2) +
  theme_cowplot() +
  xlab('max change in \n center of pressure') +
  ylab('positive\namplitude peaks')

# ggplot(sub[sub$movement_condition=='internal rotation', ], aes(x = max_COPc, y = max_amp_c_around_move))+geom_point(size = 3)+geom_smooth(method='lm', color = 'red', size = 2)+theme_cowplot()+xlab('max change in \n center of pressure')+ylab('positive\namplitude peaks')

Table S19 shows the correlation matrix for the maximum peaks in EMG (row and columns 1-4), change in center of pressure (5) and maximum vocal amplitude (6). For this table, the data only include the internal rotation condition.

Table 19: Table 20: Correlation matrix for the internal rotation condition.

	1	2	3	4	5
max. EMG pectoralis major
max. EMG infraspinatus	.11
max. EMG rectus abdominis	-.06	-.05
max. EMG erector spinae	.24**	.08	.10
max. change in center of pressure	.19*	.07	.01	.20*
max. positive amplitude expiration	.03	-.15	-.28**	-.15	.04

Table S21 shows the correlation matrix for the maximum peaks in EMG (row and columns 1-4), change in center of pressure (5) and maximum vocal amplitude (6) For this table, the data only include the external rotation condition.

Table 21: Table 22: Correlation matrix for the external rotation condition.

	1	2	3	4	5
max. EMG pectoralis major
max. EMG infraspinatus	.06
max. EMG rectus abdominis	.06	-.11
max. EMG erector spinae	.09	.04	-.01
max. change in center of pressure	.20*	.31***	-.17	.08
max. positive amplitude expiration	-.17	.22*	-.19	-.03	.05

Table S23 shows the correlation matrix for the maximum peaks in EMG (row and columns 1-4), change in center of pressure (5) and maximum vocal amplitude (6) For this table, the data only include the extension condition.

Table 23: Table 24: Correlation matrix for the extension condition.

	1	2	3	4	5
max. EMG pectoralis major
max. EMG infraspinatus	.29**
max. EMG rectus abdominis	.03	-.15
max. EMG erector spinae	.21*	.04	.23*
max. change in center of pressure	.11	.28**	-.09	.06
max. positive amplitude expiration	.00	.05	-.12	.17	.04

Table S25 shows the correlation matrix for the maximum peaks in EMG (row and columns 1-4), change in center of pressure (5) and maximum vocal amplitude (6) For this table, the data only include the flexion condition.

Table 25: Table 26: Correlation matrix for the flexion condition.

	1	2	3	4	5
max. EMG pectoralis major
max. EMG infraspinatus	.18
max. EMG rectus abdominis	-.01	.03
max. EMG erector spinae	.18	.22*	.15
max. change in center of pressure	.28**	.34***	.02	.22*
max. positive amplitude expiration	.10	.03	-.15	.07	.02

Figure 11: Correlation of posture and amplitude for internal rotation.

Figure S11 visualizes the correlation between posture (max. change in center of pressure) and expiration amplitude (positive peaks in amplitude) that was found for the voicing condition for the internal rotation condition in Table S19. In the expiration condition, however, there are no correlations between posture and expiration amplitude.

We finally explored the timing relationship of negative and positive peaks in the amplitude of the voice relative to the peak in speed of the hand. The density plot in Fig. S12 visualizes the temporal relationship between peak speed in movement and negative (left) and positive (right) peaks in the amplitude envelope, relative to the moment the moving hand reaches a peak in speed during the trial (i.e., t = 0, at peak speed of the wrist). Our findings corroborate our general observation that positive (rather than negative) peaks seem to be more robustly driven by gesture kinetics. In the timing relations obtained, the negative peaks that we observe are more variably distributed around peak speed, while for the positive peaks we see clear occurrences just around the movement onset (just before the peak speed at t = 0). This result, next to our more robust relationship between positive peaks and the muscele and posture signals, makes us conclude that in general the movements performed in our experiment are increasing subglottal pressures increasing the voice’s amplitude.

Show code producing the plot

library(RColorBrewer)
sub <- tr_wd[ (tr_wd$trialindex>9) & (tr_wd$vocal_condition=='expire'),] #only real trials, only exhalation
#make new dataset for faceting the plot
subpeak <- sub %>% 
  pivot_longer(cols = c("max_amp_time_around_move", "min_amp_time_around_move"),
               names_to = "peak_type",
               values_to = "time")

subpeak$peak_type <- recode(subpeak$peak_type,
                            max_amp_time_around_move = "positive peak",
                            min_amp_time_around_move = "negative peak")


densitydistance <- subpeak %>% 
  ggplot(aes(x = time, color = peak_type, linetype = movement_condition)) +
  geom_density(size = 0.75) +
  scale_linetype_manual(values = c("extension"="dotted", "external rotation"="dashed", "flexion"="dotdash", "internal rotation"="longdash", "no movement" = "solid")) +
  #scale_colour_manual(values = c("red", "blue"), guide = "none") +
  scale_color_brewer(palette="Dark2"
                     , guide = "none") +
  xlab("distance (in ms)") +
  ylab("density") +
  labs(linetype="movement condition") +
  theme_cowplot() +
  theme(axis.text.x = element_text(angle = 30),
        legend.key.size = unit(0.5, "cm")) +
  background_grid() +
  facet_grid(~peak_type)

Figure 12: Density plot showing temporal distance (in milliseconds) of positive and negative peaks in amplitude relative to peak speed in wrist movement.